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Abstract  of  the  Dissertation 


Particle-laden  thin  film  flow: 

An  alternating  direction  implicit  scheme  and 
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Gravity- driven  thin  film  flows  have  been  analyzed  in  terms  of  fourth-order  lubrica¬ 
tion  models,  similarity  solutions,  traveling  wave  solutions,  numerical  simulations, 
and  experiments.  However,  in  the  case  where  particle  are  suspended  within  the 
fluid,  studies  have  been  largely  limited  to  lubrication  models,  one-dimensional 
numerical  simulations,  and  experiments.  We  present  a  numerical  scheme  for  a 
lubrication  model  derived  for  particle-laden  thin  film  flow  in  two  dimensions  with 
surface  tension.  The  scheme  relics  on  an  alternating  direction  implicit  process  to 
handle  the  higher-order  terms,  and  an  iterative  procedure  to  improve  the  solution 
at  each  timestep.  Several  aspects  of  the  scheme  are  examined  for  a  test  problem, 
such  as  the  timestep,  runtime,  and  number  of  iterations.  The  results  from  the 
simulation  are  compared  to  experimental  data.  The  simulation  shows  good  qual¬ 
itative  agreement.  It  also  suggests  further  lines  of  inquiry  for  the  physical  model. 
For  constant-volume  particle-laden  thin  film  flow,  a  lubrication  model  with  pre¬ 
cursor  and  experiments  are  compared  to  a  power  law  for  the  position  of  the  front 
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of  the  flow  with  respect  to  time.  This  power-law  behavior  was  originally  derived 
for  clear  fluid  flows.  In  the  lubrication  model,  the  precursor  has  a  large  effect 
on  the  speed  of  the  front,  independent  of  the  settling  of  the  particles.  Compar¬ 
ison  between  theory  and  experiments  indicates  that  this  scaling  law  persists  to 
leading  order  for  particle-laden  thin  him  hows  with  particle  settling.  For  gravity- 
driven  particle-laden  thin  him  hows  on  an  inclined  plane,  three  distinct  regimes 
can  be  observed:  particles  settling  to  the  substrate,  a  particle-rich  ridge  forming 
at  the  front  of  the  how,  and  the  particles  staying  well-mixed.  Experiments  are 
conducted  for  a  variety  of  particle  sizes  and  liquid  viscosities.  We  compare  ex¬ 
perimental  results  with  equilibrium  theory  that  balances  shear-induced  migration 
and  hindered  settling.  We  find  that  the  well-mixed  regime  is  transient,  with  the 
particle  size  and  liquid  viscosity  influencing  its  time  scale. 
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CHAPTER  1 


Introduction 

Gravity-driven  particle-laden  thin  film  flow  occurs  in  many  different  contexts, 
from  mudslides  to  oil  spills.  A  lubrication  model  has  been  derived  for  such  flows 
in  two  dimensions  [CBH08,  CAB09].  The  model  is  formulated  in  terms  of  the 
film  thickness  and  particle  concentration,  with  the  two  dimensions  lying  in  the 
substrate  (see  Figure  2.1).  It  assumes  that  the  particles  migrate  to  the  front  of  the 
flow,  forming  a  particle-rich  ridge.  This  is  one  of  three  possible  regimes  observed 
in  experiments:  particles  settling  to  the  substrate,  migrating  to  the  front  of  the 
flow  and  forming  a  particle-rich  ridge,  or  staying  well-mixed  within  the  fluid. 
These  regimes  are  determined  by  varying  the  initial  particle  concentration  and 
inclination  angle.  Particle-laden  thin  film  flows  can  be  examined  using  numerical 
simulations  of  the  model,  experiments,  and  theory  existing  for  similar  problems 
such  as  clear  fluid  flows. 

In  recent  years,  the  problem  of  numerically  solving  gravity-driven  thin  film 
flow  for  clear  fluids  has  had  ample  work  done  in  both  one  and  two  dimensions. 
However,  the  case  when  the  film  contains  particles  suspended  within  it  has  re¬ 
ceived  less  attention,  especially  in  two  dimensions.  The  evolution  of  a  clear  fluid 
down  an  inclined  plane  is  modeled  using  a  single  partial  differential  equation  and 
numerical  schemes  have  been  derived  using  finite  differences  [DK02,  Kon03a]  and 
finite  elements  [SRX07].  For  similar  equations,  such  as  spreading  thin  films,  there 
are  methods  for  finite  elements  in  one  dimension  [GROO,  GR01,  ZBOO]  and  for 
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finite  differences  in  two  dimensions  [WB03].  The  incorporation  of  particles  into 
such  a  flow  leads  to  another  variable  in  the  model,  namely  the  particle  concen¬ 
tration,  and  an  accompanying  equation  related  to  the  evolution  of  the  particles. 
The  result  is  a  system  of  equations  that  requires  a  different  approach  from  the 
clear  fluid  case  to  formulate  a  practical  numerical  scheme,  due  to  the  coupling  of 
the  equations. 

An  active  area  of  research  in  the  last  decade  has  been  the  development  of 
numerical  methods  for  higher-order  thin  film  equations  including  complex  flu¬ 
ids  described  by  systems  of  equations.  Related  problems  include  methods  for 
coupled  systems  of  nonlinear  parabolic  equations  [LCS06,  Pao99].  The  scheme 
presented  here  is,  in  part,  inspired  by  recent  models  for  surfactants  [WCM05]  and 
thin  films  [WB03].  We  choose  an  Alternating  Direction  Implicit  (ADI)  scheme 
as  a  tractable  method  for  implicit  timesteps,  because  surface  tension  introduces 
a  severe  restriction  on  the  timestep  in  the  case  of  explicit  schemes.  This  ADI  ap¬ 
proach  also  allows  for  an  implicit  scheme  while  avoiding  to  have  to  solve  the  large 
sparse  linear  algebra  problems  by  an  iterative  method,  such  as  GMRES,  that  re¬ 
sult  from  linearizing  the  two-dimensional  operators  in  Newton’s  method  [WB03]. 
ADI  is  also  amenable  to  parallelization.  While  ADI  schemes  for  numerically  solv¬ 
ing  parabolic  equations  date  back  to  the  1950’s  [PR55],  their  use  in  higher-order 
problems  is  rather  new,  e.g.,  [WB03],  and  not  all  that  well-studied.  However, 
the  ease  of  parallelization  makes  such  schemes  a  viable  choice  for  multiprocessor 
platforms.  Since  their  inception,  ADI  schemes  have  been  extended  to  handle 
parabolic  problems  with  mixed  derivative  terms  [BW80,  CS88,  McK71,  SIJ76], 
variable  coefficients  [Kar09,  WB03],  and  high-order  terms  [WB03]. 

The  ideas  present  in  these  schemes  can  be  combined  to  create  an  efficient 
way  to  numerically  solve  the  particle-laden  thin  film  flow  equations.  The  nonlin- 
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earity  and  higher-order  terms  are  handled  in  a  similar  manner  to  Witelski  and 
Bowen  [WB03],  which  dealt  with  thin  film  equations,  and  the  remaining  terms  are 
treated  as  in  Warner  et  al.  [WCM05],  which  devised  a  semi-implicit  scheme  for 
surfactants.  This  combined  approach  is  fine-tuned  to  draw  out  better  efficiency, 
via  adaptive  timestepping  and  an  iterative  procedure  within  each  timestep.  At 
the  cost  of  the  extra  calculations  clue  to  the  iterative  nature  of  the  scheme,  the 
timestep  needed  for  stability  can  be  improved  over  recent  methods.  The  result  is 
an  efficient  method  to  simulate  the  continuum  model  in  two  dimensions. 

The  full  physics  of  particle-laden  thin  film  flow  is  not  well  understood.  Re¬ 
cent  experiments,  and  their  comparison  to  the  model,  have  raised  questions.  We 
present  such  a  comparison  in  this  paper,  where  the  results  show  qualitative  agree¬ 
ment.  In  particular,  by  performing  two-dimensional  simulations,  we  are  able  to 
observe  finger  formation  and  compare  directly  with  experiments.  The  develop¬ 
ment  of  quantitatively  correct  models  for  these  systems  is  an  ongoing  active  area 
of  research.  Thus,  there  is  a  need  for  accurate,  fully  two-dimensional  simulations 
of  the  model. 

Particle-laden  thin  film  flow,  as  treated  in  two-dimensional  simulations,  is 
most  practical  as  a  constant-flux  problem.  However,  experiments  are  done  for 
the  constant- volume  problem.  A  scaling  law  for  the  front  of  the  flow,  x(t )  Ct V3 
[Hup82]  exists  in  the  constant- volume  case  for  clear  fluids  and  our  results  indicate 
that  this  is  still  valid  for  the  particle-laden  case,  to  leading  order.  We  also  test  the 
impact  of  choosing  different  hindered  settling  functions,  which  model  the  effects 
of  particles  settling  amongst  other  particles.  For  particle-laden  flows,  we  compare 
the  results  for  theory,  numerical  simulations,  and  experiments.  This  comparison 
is  performed  using  the  scaling  constants  C  for  each  to  see  how  well  the  three 
cases  agree. 
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For  the  particle-laden  thin  him  how  lubrication  model,  we  assume  that  we  are 
in  the  regime  where  a  particle-rich  ridge  forms  at  the  front  of  the  how.  However, 
depending  on  the  initial  particle  concentration  and  inclination  angle,  different 
behaviors  can  be  observed.  We  carry  out  a  systematic  study  of  the  regimes 
for  various  particle  sizes  and  liquid  viscosities.  We  compare  experiments  with 
equilibrium  theory  to  reveal  the  transient  nature  of  the  well-mixed  regime,  where 
the  particles  neither  appear  to  settle  to  the  substrate  nor  form  a  particle-rich 
ridge.  The  experiments  also  reveal  how  the  particle  size  and  liquid  viscosity 
affect  the  time  scale  of  this  transient  regime. 

This  dissertation  is  organized  as  follows:  Chapter  2  discusses  the  model  for 
particle-laden  thin  him  how  in  two  dimensions.  In  Chapter  3,  we  present  a  full 
derivation  of  the  numerical  scheme  for  this  problem,  with  both  temporal  and 
spatial  discretizations  as  well  as  a  discussion  of  a  moving  reference  frame.  Simu¬ 
lations  and  choices  for  implementing  the  numerical  scheme  are  given  in  Chapter  4. 
The  constant-volume  problem  and  comparison  between  theory,  numerical  simula¬ 
tions,  and  experiments  are  presented  in  Chapter  5.  A  comparison  of  experiments 
to  equilibrium  theory  is  given  in  Chapter  6.  A  summary  of  the  previous  chapters 
and  results  is  presented  in  Chapter  7. 
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CHAPTER  2 


Model 

The  results  from  experiments  indicate  that  particle-laden  thin  film  flows  exhibit 
three  distinct  regimes,  based  on  the  initial  particle  concentration  and  angle  of 
inclination  [ZDB05].  For  low  concentrations  and  angles,  the  particles  settle  to  the 
substrate  with  clear  fluid  flowing  over  the  top.  The  behavior  after  sedimentation 
is  similar  to  clear  fluid  experiments,  such  as  those  performed  by  Huppert  [Hup82], 
High  concentrations  and  angles  cause  a  particle-rich  ridge  to  emerge  at  the  front 
of  the  flow.  Medium  concentrations  and  angles  lead  to  a  particle  concentration 
which  appears  to  stay  well-mixed  throughout  the  duration  of  the  experiment. 
Based  on  Cook  [C00O8] ,  this  behavior  likely  belongs  to  one  of  the  two  previously 
mentioned  regimes,  but  may  not  have  evolved  to  the  point  where  this  distinction 
can  be  made. 

The  evolution  equations  for  the  flow  are  based  on  the  regime  where  the  in¬ 
clination  angle  and  particle  concentration  are  both  high  enough  to  induce  the 
formation  of  a  particle-rich  ridge.  The  equations  are  formulated  in  terms  of  the 
thickness  of  the  film,  h,  and  the  particle  concentration  by  volume,  <f>  (see  Figure 
2.1).  The  equations  for  modeling  this  regime  were  first  derived  in  Zhou  et  al. 
[ZDB05] ;  re-derived  in  Cook  et  al.  [CBH08],  using  conservation  of  volume  rather 
than  mass;  and  modified  in  Cook  et  al.  [CAB09],  adding  in  a  shear-induced  dif¬ 
fusion  term  to  correct  for  an  instability  affecting  q h  The  dimensionless  system 
[CAB09]  is 
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Figure  2.1:  The  coordinate  system  and  variables  considered  in  this  problem,  x 
is  in  the  plane,  in  the  direction  of  the  flow;  y  is  in  the  plane,  perpendicular  to 
x;  and  z  is  normal  to  the  plane,  h  is  the  film  thickness  and  0  is  the  particle 
concentration. 


ht  +  V  •  (/ivav)  =  0,  (2.1) 

(4>h)t  +  V  ■  [<f>h  (vav  +  (1  -  0)vrel)  -  Fdiff]  =  0.  (2.2) 

The  orientation  for  (2.1)-(2.2)  is  such  that  x  lies  in  the  plane  and  is  parallel  to 
the  direction  of  the  flow,  y  is  across  the  plane  and  perpendicular  to  x,  and  z  is 
normal  to  the  plane.  A  one-dimensional  form  of  the  problem  considers  only  the 
^-direction,  while  two  dimensions  includes  both  x  and  y.  The  two  velocity  terms, 
vav  and  vrei,  are  the  volume-averaged  velocity  of  the  fluid  and  the  velocity  of  the 
particles  relative  to  the  liquid,  respectively.  We  use  the  term  liquid  to  refer  to  the 
substance  that  the  particles  are  suspended  in  and  fluid  to  refer  to  the  mixture 
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as  a  whole.  In  Equation  (2.2),  vav  +  (1  —  0)vrei  is  the  individual  velocity  of  the 
particles  [CBH08]  and  Fdiff  is  shear-induced  diffusion  of  the  particles. 

The  volume-averaged  velocity  of  the  liquid  and  the  particles  together  is 

v^m^h-D(a)  [mv{p(m-smv{pW)\ +WA’  (2'3) 

where  the  terms  in  (2.3)  come  from  surface  tension,  the  effects  of  gravity  normal 
to  the  inclined  plane,  and  the  effects  of  gravity  parallel  to  the  inclined  plane. 

The  density  of  the  fluid  as  a  whole  is  p{4>)  =  1  +  p/0;  pf  =  — — —  is  the 

Pi 

difference  in  the  densities  between  the  particles  and  the  liquid.  The  function 
p(4>)  =  (1  —  (j) / 0max)  2  [Kri72,  SP05]  is  the  effective  fluid  viscosity,  where  0max 
is  the  maximum  packing  fraction  of  particles,  assuming  the  particles  are  spheres. 
For  this  problem,  the  maximum  packing  fraction  has  been  empirically  determined 
to  be  0.58,  while  the  theoretical  value  is  0.64  [WWG09].  D(a)  =  {3Ca)1^  cot  a 
[BB97]  is  a  modified  capillary  number,  where  Ca  is  the  capillary  number  of  the 
liquid  and  a  is  the  angle  of  inclination  of  the  plane  on  which  the  fluid  is  flowing 
(a  =  0  corresponds  to  the  plane  being  horizontal  while  a  =  7t/2  to  vertical). 

The  settling  velocity  of  the  particles,  relative  to  the  velocity  of  the  liquid,  is 
a  combination  of  three  factors,  assumed  to  be  multiplicative, 

Vrel  =  Vsf((j))w(h)±  (2.4) 

2 

The  coefficient  Vs  =  -a2pf  in  (2.4)  is  the  Stokes  settling  velocity  of  a  single 

o 

sphere  settling  in  a  viscous  liquid,  where  a  is  the  dimensionless  particle  ra¬ 
dius.  A  hindered  settling  function,  in  this  case  the  Richardson- Zaki  function 
f(4>)  =  (1  —  (p)5  [RZ54a],  accounts  for  the  effect  of  sedimentation.  The  parti¬ 
cles  settling  parallel  to  the  substrate  is  modeled  using  a  wall  effects  function, 
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w{h)  =  A(h/a)2/y/ 1  +  ( A[h/a )2)2  with  A  =  1/18.  This  function  is  an  approx¬ 
imation  to  a  method  of  images  solution  to  a  single  sphere  falling  parallel  to  a 
vertical  wall  [HB65] .  This  has  the  property  that  it  is  near  0  for  h  small  and  near 
1  for  h  large. 

Since  the  system  (2.1)-(2.2)  is  fourth-order  and  (2.3)  contains  higher-order 
terms  but  (2.4)  does  not,  vrei  is  not  regularized.  This  leads  to  an  instability 
affecting  the  particle  concentration  in  numerical  simulations  [CAB09].  To  correct 
for  this,  a  shear-induced  diffusion  term  (2.5)  was  added  in, 

Fdiff  =  ^a2(3Ca)1/3I)(0)^^V0.  (2.5) 

2  Mw 

This  behavior  can  be  seen  in  a  one-dimensional  example  on  the  domain  x  :  0  — 
50  with  Ax  =  0.05.  The  initial  film  thickness  is  a  jump,  from  1  to  0.05  at  x  =  25, 
smoothed  by  hyperbolic  tangent.  The  initial  particle  concentration  is  taken  to  be 
(j)  =  0.3.  This  simulation  is  similar  to  those  described  in  Chapter  4,  and  a  moving 
reference  frame  is  used,  as  discussed  in  Section  3.4.  By  time  t  =  1000,  the  solution 
without  the  extra  diffusion  term  has  developed  an  instability  near  x  =  10  (Figure 
2.2)  while  the  one  with  it  is  still  stable  (Figure  2.3).  Note  that  the  oscillations 
trailing  the  particle-rich  ridge,  between  x  =  0  and  x  =  10,  are  a  result  of  the 
discretization  of  the  moving  reference  frame  and  are  discussed  in  Section  3.4. 
Equation  (2.5)  accounts  for  horizontal  diffusion  of  particles  in  the  fluid  caused  by 
horizontal  gradients  of  (j)  and  was  derived  based  on  results  from  Leighton  [Lei85] 
and  Leighton  and  Acrivos  [LA87b],  The  term  £)(</>)  =  (1/ 3)</>2  (l  +  (l/2)e8'8?i)  is 
a  dimensionless  diffusion  coefficient. 


Figure  2.2:  The  numerical  solution  of  (ft  at  time  t  =  1000  without  shear-induced 
diffusion.  By  this  time,  an  instability  has  developed  near  x  =  10. 


Figure  2.3:  The  numerical  solution  <ft  at  time  t  =  1000  with  shear-induced  diffu¬ 
sion  (2.5).  The  solution  is  still  stable  dne  to  the  extra  term. 
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CHAPTER  3 


Numerical  Scheme 

3.1  Time-Discretized  Scheme 

In  the  case  of  a  gravity-driven  clear  fluid  flow,  the  model  reduces  to  a  single 
equation  [BB97]  for  the  film  thickness,  h, 

ht  +  ( h3)x  +  V  •  (h3VV2/i  -  D(a)h3Vh )  =  0.  (3.1) 

Solving  (3.1),  and  similar  problems,  numerically  in  one  and  two  dimensions  has 
been  performed  using  several  different  methods  [BBG98,  DK02,  Kon03a,  LGP02, 
SRX07,  WB03].  Including  particles  in  the  physics  not  only  adds  a  second  equa¬ 
tion,  but  couples  it  to  the  equation  for  the  film  thickness.  The  particle-laden  case 
has  been  solved  numerically  in  one  dimension  with  methods  such  as  forward  Euler 
with  upwind  differencing  [ZDB05]  and  the  Lax-Friedrichs  method  [CBH08]  when 
the  high-order  terms  are  omitted,  and  backward  Euler  with  centered  differencing 
[ZDB05]  when  the  terms  are  included. 

This  system  of  PDEs  in  two  dimensions  poses  numerical  difficulties  beyond 
those  present  in  the  clear  fluid  problem.  For  both  the  clear  and  particle-laden 
cases,  fully  explicit  schemes  typically  have  the  problem  that  an  0( Ax4)  timestep, 
assuming  Arc  =  Ay,  is  needed  for  stability.  One  solution  is  to  use  an  implicit 
scheme.  For  the  clear  fluid  and  similar  problems,  the  nonlinearity  combined  with 
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an  implicit  scheme  amounts  to  solving  the  problem  at  each  timestep  using  an 
iterative  process,  such  as  Newton’s  method,  to  converge  to  the  solution  [WB03]. 
For  the  particle-laden  case,  using  an  implicit  scheme  typically  requires  that  both 
equations  be  solved  simultaneously,  using  an  iterative  process  to  account  for  the 
nonlinearity.  This  results  in  a  linear  algebra  problem  with  twice  the  number  of 
unknowns  and  a  matrix  that  is  twice  as  large  in  each  dimension,  compared  to 
the  clear  fluid  problem.  Therefore,  solving  the  particle-laden  case  leads  to  larger 
linear  algebra  problems  to  solve  at  each  timestep  and  the  matrix  from  Newton’s 
method  will  have  a  more  complex  structure  than  for  clear  fluids. 

The  goal  of  the  scheme  presented  here  is  to  circumvent  some  of  the  afore¬ 
mentioned  difficulties.  The  advantages  of  this  approach,  over  a  purely  explicit 
scheme  or  implicit  with  Newton’s  method,  are  that  the  timestep  is  more  lenient 
than  for  a  fully  explicit  scheme  and  the  linear  algebra  problem  that  results  from 
the  implicit  part  of  the  scheme  is  reduced  to  a  series  of  smaller  banded  matrix 
solves,  which  can  be  done  efficiently  and  independently  for  each  equation. 

The  numerical  scheme  that  we  employ  for  the  particle-laden  thin  film  flow 
problem  is  inspired  by  the  schemes  presented  in  Witelski  and  Bowen  [WB03]  for 
higher-order  parabolic  PDEs  and  Warner  et  al.  [WCM05]  for  surfactants.  In 
Witelski  and  Bowen,  several  ADI  schemes,  based  on  backward  Euler,  second- 
order  backward  difference  formulas,  as  well  as  Newton-like  schemes,  are  derived 
for  solving  the  nonlinear  PDE  known  as  the  thin  film  equation, 

ht  +  V-  (/(/j)VV2/j)  =  0.  (3.2) 

The  backward  Euler-based  ADI  scheme  for  (3.2)  uses  approximate  values  of  h  in 
the  nonlinear  and  mixed-derivative  implicit  terms.  It  is  suggested  to  start  with 
approximations,  such  as  time-lagged  values,  for  evaluating  these  terms  and  calcu¬ 


li 


lating  the  numerical  solution  at  the  timestep.  Then  use  this  solution  for  the  new 
approximate  values  within  the  same  timestep  and  recalculate.  This  results  in  an 
iterative  scheme  at  each  timestep.  However,  for  solving  the  thin  him  equation,  it 
was  noted  that  the  iterations  did  not  provide  a  noticeable  improvement.  Warner 
et  al.  use  this  method  for  a  coupled  system  of  nonlinear  PDEs  relating  to  sur¬ 
factants.  They  handle  the  higher-order  terms  implicitly  using  Crank- Nicolson, 
and  apply  ADI  to  this.  The  remaining  terms,  which  are  at  least  second-order 
in  space,  are  treated  explicitly.  For  the  nonlinearity  and  mixed-derivative  terms, 
the  values  are  time-lagged  and  the  problem  is  solved  only  once  per  timestep.  In 
the  simulations,  Ax  =  Ay  =  7r/100  ~  0.0314  required  a  timestep  of  0(1CT5). 

Our  approach  is  to  handle  applicable  terms  implicitly,  using  ADI,  and  treat 
the  remaining  terms  explicitly,  as  we  show  below.  The  terms  handled  implicitly 
are  those  with  spatial  derivatives  on  the  same  variable  as  the  time  derivative. 
For  example,  Equation  (2.1)  has  the  time  derivative  on  h,  so  the  terms  treated 
implicitly  should  have  spatial  derivatives  on  h.  Making  this  choice  allows  for  the 
splitting  of  the  two-dimensional  operators  into  the  product  of  two  one-dimensional 
operators  in  the  derivation  of  the  ADI  scheme.  Iterations  within  each  timestep 
allow  for  a  larger  At  to  be  taken  at  the  cost  of  some  extra  calculations.  In  general, 
the  increase  in  the  size  of  the  timestep  outweighs  the  extra  computational  work, 
as  shown  in  Chapter  4. 

For  Equation  (2.1),  the  terms 


V- 


VV2h  + 


(3.3) 


can  be  handled  implicitly.  This  is  because  the  spatial  derivatives  on  these  terms 
are  applied  to  h.  Of  these  terms,  some  parts  of  them  will  be  handled  by  ap¬ 
proximation,  as  in  Witelski  and  Bowen  [WB03].  Including  the  first-order  terms 
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in  the  implicit  treatment  allows  them  to  be  discretized  spatially  using  centered 
differencing  to  maintain  stability.  Solving  this  equation  numerically  assumes  that 
(j)  is  known,  or  can  be  approximated,  and  we  are  solving  for  h.  First  discretize 
the  terms  in  (3.3)  in  time  with  backward  Euler,  including  the  time  derivative, 

hn+1  +  AtV  ■  ( -^-VV2h+ ^-h3±]  =hn.  (3.4) 

VM0)  Mw  J 

Write  out  the  operators  in  (3.4)  fully, 

hn+1  +  At  dx  (  .  hxxx]  +  9y  (  hyyy\ 

[  Vh(0)  J  Vh(0)  J 

+dx  +  At  dx  (w)hyyx)  +  dy  (w)hxxy).  = 

The  idea  behind  the  ADI  approach  is  to  reduce  the  implicit  part  of  (3.5), 
with  derivatives  in  both  x  and  y,  to  a  product  of  two  operators,  each  with  only 
derivatives  in  either  x  or  y.  To  achieve  this,  the  terms  involving  only  ^-derivatives 
and  only  ^-derivatives  are  grouped  together.  Define  the  operators 

Dx  =  dx  (x^)dxxx  +  W)h21)  r  Dy  =  dy  {jmdyyy)  '  (3-6) 

Then  replacing  the  terms  in  (3.5)  with  the  definitions  in  (3.6),  we  have 

hn+1  +  A  t(Dx  +  Dy)hn+1 

r  /  h3  \  (  h3  \in+1 

+At  vx  \Mhyyx)  +  dy  \lM)hxxv ) \  =  }>n'  (3'7) 

In  order  to  obtain  an  ADI  scheme  from  (3.7),  note  that  /  +  A tDx  +  A tDy  = 
(/  +  AtDx)(I  +  AtDy)  —  (A t)2DxDy  and  so  the  left-hand  side,  with  the  addition  of 
an  O  (At2)  term,  can  be  written  as  a  product  of  two  one- dimensional  operators. 
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(/  +  A tDx)(I  +  A tDy)hn+1  -  (A tfDxDyhn+1 

r  /  h3  \  (  h3  \in+1 

+At  [dx  \jM)hyyx)  +  dy  \lM)hxxy ) \  =  h'n'  (3'8) 

To  handle  the  nonlinear  terms,  which  occur  in  front  of  derivatives,  and  mixed- 
derivative  terms  in  (3.8),  define  them  as  approximate,  denoted  by  a  tilde  (e.g., 
hn+l).  The  approximate  terms  can  be  chosen  in  some  reasonable  manner,  such  as 
time-lagged  or  extrapolated.  This  will  be  discussed  in  more  detail  later.  Subtract 
the  mixed-derivative  terms  from  and  add  the  0(Af2)  term  to  both  sides.  This 
leaves  a  scheme  in  which  all  the  terms  operating  on  hn+1  are  known,  as  is  the 
entire  right-hand  side. 

(/  +  A  tDx)(I  +  A  tby)hn+l  =  hn 

+  l^(At)2DxDy  -  At  dx  +  dy  |  hn+1.  (3.9) 

For  simplicity,  define  the  operators  in  (3.9)  as 

Lx  =  I  +  A tDx,  Ly  =  I  +  A tDy. 

Subtracting  LxLyhn+l  from  both  sides  of  (3.9),  which  cancels  the  0(At2)  term, 
yields 

LxLy  (V+1  -  /A+1)  =  -  (hn+1  -  hn^j 

~  \  n+1 

-tw2H%x  .  (3.10) 

v(4>)  d(4>)  J 
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At  this  point,  the  implicit  part  of  the  scheme  is  complete  and  the  explicit  terms 
can  be  added  back  into  (3.10)  using  forward  Euler. 


+AfV  •  <  D(a 


LxLy  (hn+1  -  hn+1^j  =  -  (hn+1  -  hn^j 

~  \  n+1 

^tvv2n^ft3x 


(3.11) 


Define 


u  =  hn+1  -  hn+l , 

which  can  be  thought  of  as  a  correction  term  to  the  approximation  of  hn+1,  and 
(3.11)  can  be  written  as  a  three-step  process:  two  one-directional  solves  (3.12)- 
(3.13)  and  an  update  step  (3.14). 


Lxv  =  —  { hn+l  —  hn 


A tV  ■ 


~h  -VV2H^t/i3x 


n+1 


+AtV  •  <  D(a 


(3.12) 


LyU  =  V, 


(3.13) 


/i”+1  ss  hn+1  +  u.  (3.14) 

Since  the  operators  Lx  and  Ly  involve  at  most  fourth-order  terms,  the  spatial 
discretization  of  them  will  lead  to  a  five-point  stencil  in  the  x-  and  ^/-direction, 


15 


respectively.  This  discretization  is  discussed  fully  in  Section  3.3.  Along  each 
row/column  of  the  discretized  domain,  this  results  in  a  pentadiagonal  linear  alge¬ 
bra  problem.  This  can  be  solved  using  a  pentadiagonal  solver,  or  a  more  generic 
banded  matrix  solver. 

To  help  with  the  inaccuracy  in  the  nonlinear  and  mixed-derivative  terms  re¬ 
sulting  from  approximation,  an  iterative  procedure  can  be  used  at  each  timestep 
to  improve  the  solution  and  size  of  the  timestep.  This  was  first  suggested  for 
the  ADI  scheme  in  the  context  of  thin  film  equations  [WB03].  This  procedure 
amounts  to  repeating  the  three-step  process  associated  with  solving  each  equa¬ 
tion  at  each  timestep  and  updating  the  approximate  solution  with  the  most  recent 
solution,  until  the  new  and  approximate  solutions  sufficiently  converge.  For  ex¬ 
ample,  one  would  solve  (3.12)-(3.14),  solve  (3.24)-(3.26),  and  examine  how  much 
the  approximate  solution  differs  from  this  computed  solution.  If  this  difference 
is  significant,  one  can  replace  the  old  approximate  terms  with  the  computed  so¬ 
lution  and  solve  the  same  timestep  again.  This  process  can  be  continued  until 
the  approximate  and  computed  solutions  are  close.  This  is  similar  to  fixed-point 
iteration. 

For  Equation  (2.1),  when  entering  the  timestep,  a  choice  must  be  made  as  to 
the  value  of  hn+ 1  and  (<ph)n+1.  Using  h  as  an  example,  two  reasonable  choices 
would  be  a  time-lagged  approximation,  hn,  which  is  a  first-order  accurate  ap¬ 
proximation  in  time,  or  an  extrapolated  approximation,  2 hn  —  hn_1,  which  is 
second-order  in  time.  For  adaptive  timestepping,  this  extrapolation  is  given  by 
hn  +  (At/  At  0\d)(hn  —  hn_1),  where  At  is  the  prospective  timestep  between  tn  and 
tn+ 1  and  At0id  is  the  timestep  between  tn_1  and  tn.  While  the  second  choice  of 
an  approximation  is  second-order,  it  also  requires  storing  an  extra  set  of  data, 
namely  hn~l.  Other  choices  for  estimating  hn+1  and  (<ph)n+1  based  on  previous 


16 


data  could  be  used  as  well.  With  this  choice  made,  the  three-step  process  for 
each  equation  can  be  implemented,  obtaining  a  solution  for  hn+1  and  (0h)n+1. 
We  refer  to  the  case  when  the  solution  obtained  here  is  accepted  as  perform¬ 
ing  One  Iteration.  However,  at  this  point,  the  approximation  can  be  redefined, 
hn+1  =  hn+1,  and  the  process  repeated.  This  can  be  continued  until  convergence 
between  the  approximate  and  new  solution,  or  equivalently  when  the  correction 
term  u  is  small  in  a  chosen  norm.  We  refer  to  this  case  as  Iteratioiis  since  the 
problem  is  solved  iteratively  for  each  timestep. 

For  (2.2),  the  ADI  method  is  applied  to  (j)h  as  a  whole,  since  the  time  derivative 
is  on  this  term.  The  applicable  terms  in  the  equation  are 


+4>h 


P(0) 

M0) 


ft2  +  (1  -  <p)V,f(<t>)w(h) 


(3.15) 


As  with  (2.1),  the  time  discretization  of  (3.15)  is  based  on  a  backward  Euler 
method 


+(J)h 


(0h)n+1  +  AfV  • 
P(<f>) 


-D(a)  I  V(0h) 


P{<t>) 


p{4>) 


h 2  +  (1  -  4>)Vaf(<f>)w{h) 


1  n+ 1 


=  (#)' 


(3.16) 


Writing  out  the  operators  in  (3.16)  explicitly, 


(0h)n+1  -  A tD(a)pf 


d *  ^mdMh)) +  d> 


+A  tdx 


<ph 


p(4>) 

P{<t>) 


h2  +  (1  -  (j))Vsf((j))w(h) 


71-1-1 


=  (0h)n.  (3.17) 
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Define  the  operators  in  (3.17)  involving  only  ^-derivatives  and  only  ^-derivatives 
as  Vx  and  Vy,  respectively. 


^  ^  a  ( ((j)h)h2  a 

Vx  -  D(a)pfdx  ^  ft 


n+ 1 


+ft 


p{4>) 


h2(l  -  <t>)Vsf{<j>)w{h ) 


Vy  =  -D{a)pfd% 


( (4>h)h2 


ft 


I 

n+ 1 


71+1 


y  V  y 

Using  (3.18),  the  equation  can  be  compactly  written  as 


(3.18) 


(0h)n+1  +  At  (: Vx  +  Vy)  (0h)n+1  =  (<f>h)n.  (3.19) 

Note  that  there  are  no  mixed- derivative  terms  to  handle  in  (3.19).  The  left-hand 
side  can  be  written  as  the  product  of  two  one-dimensional  operators,  incurring 
an  0(At2)  term  in  the  process. 


(/  +  A tvx)  (/  +  A tvy)  (<f>h)n+1  -  (A t)2VxVy((j>h)n+1  =  (0h)n.  (3.20) 

Add  the  0(Af2)  term  to  both  sides  of  (3.20),  and  make  all  terms  that  occur 
nonlincarly  at  time  tn+l  approximate,  as  before. 


(i  +  A tV^j  (i  +  A tVy'j  (4>h)n+1  =  (4>h)n  +  (A t)2VxVy(<t>h)n+1.  (3.21) 


Define 


Cx  =  I  +  A tVx,  ty  =  I  +  A tVy 
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and  subtract  Cx£y((f>h)n+1  from  both  sides  of  (3.21)  to  obtain 


Cx£y  ((#)"+1  -  (0h)"  +  1)  =  -  (( 4>h)n+1  -  (0 hy 

-A tV  ■ 


■■D(«)  (  Pf— 

M0)  / 

+#•  (  f>^rh2  +  (1  -  <f>)Vsf{4>)w{h)  l  X 
\PW) 


-i  n+l 


(3.22) 


The  remaining  terms  can  be  incorporated  into  (3.22)  via  forward  Euler. 


CXty  (( <j)h)n+1  ~  (0h)"  +  1)  =  -  ((0h)n+1  -  (#)”) 

-A tv  • 


+<fih 


P{$) 

U(0) 


-D(a)  | 

p(4>) 

h2  +  (1  -  4>)Vsf(4>)w(h) 


n+l 


—A  tV- 


h 2 

4>h  (  — VV2h  -  £>(a) 


h(0) 


/  k2  _ ,  5  h3  _  .  .  , . . 

I  - Vh  -  (p(0)) 


V/l(0) 


8/i(0) 


1  diff 


(3.23) 


Dehne 


w  =  (0h)n+1  —  (0h)n+1. 

Then  (3.23)  can  be  written  out  as  the  three-step  process  (3.24)-(3.26): 


19 


Cxv  =  -  (( 4>h)n+1  -  (<j>h) 

-A*v-  -D(a)  Lf^^v(i'h)\ 

+4>h  (  +  (1  -  4>)Vsf(4>)w(h)  |  X 

\M0) 


n+1 


-AfV  • 


<ph 


h 2 

M0) 


vv2/i-d(«; 


/  h2  5  h3 

•  - V/i  -  (p(0)) 


V/i(0) 


8/i(0) 


—  F  Hi 


diff 


(3.24) 


=  v,  (3.25) 

(0/i)n+1  «  (#)n+1  +  w.  (3.26) 

The  spatial  operators  in  the  and  Cy  terms  are  at  most  second-order,  and 
spatial  discretization  leads  to  a  three-point  stencil  in  each  direction.  Similar  to 
(3.12)  and  (3.13),  a  tridiagonal  solver  or  banded  matrix  solver  can  be  used  to 
solve  along  each  row/column. 

Solving  the  system,  as  a  whole,  at  each  timestep  can  be  then  achieved  by 
solving  (2.1)  using  (3.12)-(3.14)  for  hn+1,  solving  (2.2)  using  (3.24)-(3.26)  for 
(0/i)n+1,  then  recovering  the  particle  concentration  as  (j)n+1  =  {<ph)n+1  / hn+1 .  Note 
that  each  solve  only  uses  values  hn ,  hn+1 ,  (f>n ,  and  0n+1,  all  of  which  are  known. 
This  scheme  can  be  solved  in  other  possible  ways.  One  might  choose  to  use, 
after  solving  (2.1),  hn+1  in  lieu  of  an  approximation  for  hn+l  for  solving  (2.2). 
Alternatively,  the  equations  could  be  solved  in  the  opposite  order. 
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3.2  Adaptive  Timestepping 


We  use  an  adaptive  timestepping  scheme  to  advance  the  solution.  The  scheme 
utilizes  the  solution  at  consecutive  timesteps  tn_1,  tn,  tn+l .  Based  on  a  measure  of 
error,  it  decides  whether  or  not  to  accept  the  new  solution,  and  if  it  is  reasonable 
to  increase  the  size  of  the  timestep.  This  is  a  modification  of  the  scheme  used  in 
Bertozzi  et  al.  [BBD94],  in  which  it  serves  as  an  estimate  of  a  dimensionless  local 
truncation  error  in  time.  Consider  the  solution  of  the  film  thickness,  h,  at  times 
and  tn+1.  Calculate  en+1  =  ( hn+1  -  hn)/hn  and  en  =  ( hn  -  hn~x)lhn. 
The  modification  from  the  original  method  is  to  divide  by  the  value  hn  at  each 
point  rather  than  =  max,; ^ { hvi:j } ,  since  it  produces  abetter-working  adaptive 

scheme  for  this  problem.  Denote  the  timestep  going  from  time  tn  to  tn+1  as  At 
and  from  tn~x  to  tn  as  At0^.  Then  define 


Error 


e 


n+1 


At 

At0\d 


(3.27) 


This  provides  a  dimensionless  estimate  of  the  local  truncation  error  in  time, 
accumulated  over  the  grid.  The  solution  will  be  accepted  if  this  error  is  less 
than  some  tolerance,  denoted  Tob.  If  the  error  is  less  than  a  smaller  tolerance, 
Tol2  <  Tob,  for  a  fixed  number  of  steps,  the  timestep  is  increased  by  a  scale 
factor.  If  the  error  is  larger  than  Toll,  the  maximum  number  of  iterations  within 
a  timestep  is  surpassed,  or  the  solution  becomes  negative,  the  timestep  is  reduced 
by  a  factor  of  2.  An  example  for  Toll  and  Tol2  would  be  10-'  x  (  Area  of  Domain) 
and  1CT9  x  (  Area  of  Domain)  respectively,  where  the  difference  in  the  tolerances 
is  at  least  an  order  of  magnitude  apart  to  prevent  the  error  from  alternating 
between  too  large  to  accept  and  small  enough  to  increase  the  timestep.  The  form 
of  these  tolerances  were  chosen  to  make  it  convenient  for  various  size  domains 
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without  having  to  change  the  tolerances  manually  for  each  domain. 

Since  (3.27)  only  takes  into  account  one  of  the  two  variables,  this  error  can 
be  computed  for  0h,  or  merely  0,  as  well.  These  two  errors  can  be  combined  into 
an  overall  measure  of  the  error  by  taking  the  maximum  of  the  two,  or  by  some 
other  reasonable  combination  such  as  adding  the  two  errors  together  or  choosing 
a  separate  set  of  tolerances  for  each. 

3.3  Spatial  Discretization 

We  use  centered  finite  differences  for  all  spatial  discretizations.  Using  the  nota¬ 
tion,  hi+i/2j  ~  (hij  +  hi+ ij)/2,  the  fourth-order  term  in  (2.1)  is 


v'(iSvv2,t 


hJ 


hi 


--WAj  ,  _  hi~l/2,  j  ,  \  /A 

^(0J+1/2J)  Mwj)  W-1/27  7 

^+1/2, j  h  _  h  ^  /At 

^(0m/2j)  W*>i+1/2J  htti- i/2j)  raa:^1/2’V  7 


hi,j+ 1/2  j  hi,j- 1/2  ,  \ 

~TL  Uhra:i/,i,j+l/2  “TT  7^xxy,i,j- 1/2  I 

u(0ip+i/2)  KVij- 1/2)  y 

,  ^l,j+ 1/2  ,  ^i,j- 1/2  u  ^  /Ani 

7  l  7777a  \llyyy,i,j+ 1/2  “TT  \llyyy,i,j- 1/2  1  /^2/- 

VM0M+1/2)  hVPij- 1/2)  y 


(3.28) 


Here,  the  third  derivatives  are  calculated  at  half-grid  points  by  differencing  con¬ 
secutive  standard  second-order  approximations.  Two  representative  examples 
are 
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h/xxx,i-\-l/2,j  ~  ( /p+2,j  “1“  / Ax  , 

hxxy,i,j +1/2  ~  ((^i+lj+1  T  /lj— lj+l)/Ax 

~(hi+ i,j  -  2 hij  +  1}j)/Ax2) /Ay. 


(3.29) 


The  two  second-order  terms  are  discretized  as 


l3 

^+1/2,  j 


^(01+1/2,. 


■  (5) vw)A).  M 


^-1/2, j 
p{(t)i-l/2,j 


/Aa:2 


z.3 

^J+l/2 

P(4>ij+ 1/2, 


/i? 


hi-1/2 


PifiiJ- 1/2, 


■( P(4>i,j)hi,j  ~  p{<Pi,j-i)hij-i)  /Ay2, 


(3.30) 


v,wVW)AJ 


7?4 

^+l/2,i 

P(4>i+l/2,j 


~(p((fii+ l,j)  P(0i,j)) 


/?4 

__(Vl/2J_ 

p{(t)i-l/2,j 


(P(<l>i,j)  ~  P(<t>i~l,j))  /Ax2 


hA 

ni,j+ 1/2 
M0iJ+l/2, 


■(p(0i,j+l)  P(0*  J  )  ) 


/if 


id- 1/2 


hi.4)i,j-l/2. 


-P(0id-l))  AV- 


(3.31) 
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The  advective  term  is  discretized  using  a  standard  centered-differencing  scheme. 

The  terms  in  (2.2)  are  discretized  in  the  same  manner  since  many  of  them  are 
similar  to  those  in  (2.1).  The  fourth-  and  second-order  terms  that  come  from  vav 
are  discretized  as  in  (3.28)-(3.31),  with  h  replaced  by  (j)h.  Both  advective  terms 
are  discretized  via  standard  centered  differencing. 

The  shear-induced  diffusion  term  is  discretized  the  same  way  as  (3.30)-(3.31). 


V  • 


h2p((p) 

/i(0) 


V</> 


rvj.  ^i+l/2jP(^+l/2j)  (M  x  x 

D\(j)i+ 1/2 j)  "( x  \  vPi+ lj  Vhj) 


*+1/2, 


A/,  x  h?_1/2jP(<f>i-i/2,j)  .  2 

1/2, j)  77  7  (0ij  0i— lj)  ) 

p{(Pi-l/2,j) 


p{4>i 


i, 3+1/2) 


A/,  'h'ij-1/2p((f)i,j-l/2)  .  2 

1/2)  n^cf)-  y~)  (^*4  1)  )  7^2/ 


(3.32) 


Centered  differencing  is  not  used  for  the  moving  reference  frame,  if  one  is 
employed.  Instead,  a  second-order  upwind  differencing  scheme  is  used,  which 
will  be  discussed  in  the  next  section. 


3.4  Reference  Frame 

The  area  of  interest  in  the  simulations  is  near  the  front  of  the  flow,  where  effects 
like  the  capillary  and  particle-rich  ridges  occur.  With  a  fixed  reference  frame, 
the  spatial  domain  would  need  to  be  taken  as  the  entire  area  over  which  the 
flow  would  evolve,  leading  to  large  portions  of  the  domain  where  no  change  is 
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occurring.  This  issue  can  be  easily  addressed  by  using  a  moving  reference  frame. 

To  implement  a  moving  reference  frame,  we  add  an  extra  term  to  each  equa¬ 
tion,  — shx  on  the  left-hand  side  of  (2.1)  and  —s((f>h)x  on  (2.2).  ffere,  s  >  0 
is  the  speed  at  which  the  moving  reference  frame  travels.  Zhou  et  al.  [ZDB05] 
approximate  the  front  speed  by  removing  all  terms  from  the  equations  which  are 
higher  than  first  order,  leaving  only  the  advective  terms.  They  observe  that  these 
terms  capture  the  large  scale  dynamics,  including  the  speed  of  the  shocks,  and 
the  ridges  that  develop  in  h  and  <fi.  This  leaves  a  2  x  2  system  of  conservation 
laws  of  the  form 


ht  +  \F(h,  4>h)]x  —  0,  (3.33) 

(<f>h)t  +  [G(h,(j>h)]x  =  0,  (3.34) 

where  F  and  G  are  defined  as 


M0) 


p{4>) 


G(h,  4>h)  =  —-FUh)h2  +  (<j>h)(  1  -  <j>)Vsf(<j>)w(h). 
PVP) 


The  initial  conditions  for  (3.33)-(3.34)  are 


h(x,  0) 


hi,  x  <  0, 
hr,  x  >  0, 


(3.35) 


{4>h){x,  0) 


<!>ohh  x  <  0, 
0o hr,  x  >  0. 


(3.36) 


25 


where  hi  and  hr  in  (3.35)  and  (3.36)  are  the  initial  him  thickness  and  the  height  of 
the  precursor  b ,  respectively,  and  0o  hi  (3.36)  is  the  initial  particle  concentration 
of  the  fluid.  These  initial  conditions  specify  a  Riemann  problem  [Lax73].  From 
the  initial  shock  in  both  equations,  an  intermediate  state  emerges,  (hi,  (<ph)i). 
The  weak  form  of  this  system  produces  two  Rankine-Hugoniot  jump  conditions, 
which  define  the  shock  speeds,  ahead  and  behind  the  intermediate  states.  For  si, 
the  speed  of  the  shock  behind  the  intermediate  state,  and  s2,  the  speed  ahead, 
these  conditions  are  given  by 


F(hj,  (#0,)  -  F(hh  (#0;)  _  G(hj,  (0h)t)  -  G(hh  (0fr)z) 
hi  —  hi  (4>h)i  —  (4>h)i 

F(hr,  (<j>h)r)  -  F(hj,  (<j>h)j)  =  G(hr,  (<f>h)r)  -  G(hj,  (<j>h)i) 

hr  ~  hi  ( (f>h)r  -  (4>h)i 


Figure  3.1:  The  intermediate  states  that  develop  in  the  him  thickness  (left)  and 
particle  concentration  (right)  for  the  first-order  system  of  equations. 

The  intermediate  states  and  shocks  can  be  seen  in  Figure  3.1.  This  nonlinear 
system  (3.37)  of  four  equations  and  four  unknowns,  hi,((j)h)i,  s\,  and  S2,  can 
be  solved  via  Newton’s  method.  For  the  simulations  shown  in  Chapter  4,  our 
reference  frame  speed  is  an  average  of  the  two  speeds,  s  —  (si  +  s2)/ 2. 
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The  discretization  of  the  terms  for  the  moving  reference  frame  is  done  explic¬ 
itly  using  forward  Euler  combined  with  second-order  upwind-differencing, 

7  ^  „  ^h+2  ,j  T  4/ij+ij  3 

Slbnr  5  "  . 

2  Ax 

This  was  chosen  over  explicit  first-order  upwind  and  implicit  centered  differencing. 
For  a  test  run  to  time  t  —  10  with  no  variation  in  the  ^/-direction,  implicit 
centered  differencing  produced  the  highest  particle-rich  ridge,  but  introduced 
small  oscillations  ahead  of  the  flow  that  were  approximately  2%  of  the  height  of 
the  ridge.  First-order  upwind  was  dissipative  and  lead  to  the  ridge  being  28% 
smaller  than  implicit  centered  differencing.  The  effects  of  choosing  second-order 
upwind  appear  to  be  some  minor  dissipation,  about  17%  as  compared  to  implicit 
centered  differencing,  and  dispersion,  which  was  not  seen  in  this  test  case,  behind 
the  particle-rich  ridge. 

The  moving  reference  frame  can  be  used  for  both  the  one-  and  two-dimensional 
cases  (see  Figures  3.2  and  3.3).  To  demonstrate  this,  simulations  were  run  un¬ 
der  the  same  conditions  as  those  in  Chapter  4.  The  theory-based  solution  for 
the  problem  without  higher-order  terms  (3.33)-(3.36)  aligns  well  with  the  one¬ 
dimensional  numerical  solution  for  the  full  problem.  The  two-dimensional  solu¬ 
tion  for  the  full  problem  with  a  perturbation  to  the  initial  film  thickness  leads  to 
a  finger  that  moves  faster  than  the  one- dimensional  case  and  the  troughs,  to  the 
sides  of  the  finger,  move  slower. 

This  can  be  viewed  more  succinctly  in  Figures  3.4  and  3.7,  where  the  contours 
of  the  perturbed  two-dimensional  case  are  shown.  The  position  of  the  finger  runs 
ahead  of  the  one-dimensional  front,  which  is  approximately  at  x  =  15,  while  the 
troughs  lag  behind.  The  averaging  of  the  front  position  was  first  investigated  by 
Huppert  [Hup82]  for  experiments  involving  clear  fluids.  Both  simulations  start 
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with  the  same  volume  and,  after  an  initial  transient,  the  average  front  positions 
for  the  one-  and  perturbed  two-dimensional  case  (measured  at  h  =  0.5)  stay  close 
to  each  other  (Figure  3.5).  Figure  3.6  shows  the  position  of  the  finger  and  the 
trough  in  the  two-dimensional  case  over  time. 
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Figure  3.2:  Comparison  of  theory  and  simulations  at  time  t  =  100  for  the  him 
thickness,  /?,:  theory  without  higher-order  terms  (solid  line),  one-dimensional  so¬ 
lution  to  the  full  problem  (dashed  line),  perturbed  two-dimensional  finger  (dotted 
line),  and  perturbed  two-dimensional  trough  (dot-dashed  line).  The  domain  in 
the  ^/-direction  is  15  units  long,  with  the  finger  slice  taken  at  y  =  7.5  and  the 
trough  slice  taken  at  y  —  1  (see  Figure  3.4). 
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Figure  3.3:  Comparison  of  theory  and  simulations  at  time  t  =  100  for  the  particle 
concentration,  0.  The  labels  are  the  same  as  in  Figure  3.2. 


Figure  3.4:  A  contour  plot  of  the  simulation  at  times  t  =  0  (left)  and  t  =  100 
(right)  for  the  him  thickness,  h,  in  the  perturbed  two-dimensional  case.  The 
perturbation  in  two  dimensions  leads  to  a  fingering  instability  not  seen  in  the 
one-dimensional  case. 
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Figure  3.5:  The  average  front  position  of  the  him  thickness,  h,  of  the  one-dimen¬ 
sional  and  perturbed  two-dimensional  case  up  to  time  t  =  100.  After  an  initial 
transient,  the  average  front  positions  stay  close  to  each  other. 


Figure  3.6:  The  front  position  of  the  him  thickness,  h,  of  the  perturbed  two-di¬ 
mensional  case  up  to  time  t  =  100  along  the  finger  and  trough. 
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Figure  3.7:  A  contour  plot  of  the  simulation  data  at  time  t  =  100  for  the  particle 
concentration,  (p,  in  the  perturbed  two-dimensional  case.  The  perturbation  leads 
to  a  particle-rich  ridge  that  outlines  and  begins  to  fill  in  the  finger. 
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CHAPTER  4 


Simulations  and  Comparison  to  Experiments 

4.1  Benchmark  Simulations 

A  rectangular  domain  is  used  with  the  x-direction  oriented  down  the  inclined 
plane  and  the  ^/-direction  across  the  inclined  plane.  In  all  cases,  the  particle 
concentration  is  initially  taken  to  be  <f>(x,y,  0)  =  d o,  where  0  <  do  <  dmax-  This 
corresponds  to  having  a  well-mixed  initial  fluid.  The  film  thickness  far  behind 
the  contact  line  is  set  at  h(x,  y,  0)  =  1  and  ahead  of  the  flow,  a  precursor  of  height 
h(x,  y,  0)  =  b  is  assumed.  At  the  contact  line,  a  perturbation  to  a  linear  front  can 
be  applied  to  induce  behavior  such  as  a  fingering  instability.  The  parameters  in 
the  model  are  taken  to  be:  a  =  0.1,  pf  =  1.7,  Ca  =  10~3,a  —  tt/4.  The  constant 
0m ax  is  taken  to  be  0.67,  in  line  with  the  simulations  in  Cook  et  al.  [CAB09]. 
The  initial  timestep  is  set  to  At  =  10-6  and  the  mesh  width  is  Ax  =  Ay  =  0.05. 

For  the  model,  two  sources  contribute  to  the  height  of  the  film  thickness  and 
particle  concentration  near  the  front  of  the  flow.  The  first  is  the  higher-order 
terms,  such  as  surface  tension,  which  produce  smooth  ridges  in  both  h  and  <fi. 
Second,  even  without  these  terms,  an  intermediate  state  at  the  front  emerges  for 
both  variables,  higher  than  either  of  their  respective  left  or  right  states.  These 
heights  are  dependent  on  the  precursor  b. 

The  height  of  the  precursor  in  the  following  simulations  is  chosen  to  be  the 
same  as  Ax.  In  general,  the  choice  of  precursor  has  a  small  effect  on  the  speed  of 
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the  flow,  but  a  large  effect  on  both  the  film  thickness  and  particle  concentration. 
To  illustrate  this,  Table  4.1  shows  the  height  of  the  intermediate  states  for  both 
h  and  0  as  well  as  the  speeds  of  the  trailing  and  leading  shocks  obtained  from 
the  theory-based  solution  to  the  system  of  conservation  laws  (3.33)-(3.36)  (see 
Section  3.4  for  a  more  in-depth  discussion). 


b 

hi 

4>i 

Sl 

s2 

0.1 

1.01653 

0.307566 

0.459323 

0.510221 

0.05 

1.03478 

0.315538 

0.459314 

0.483782 

0.025 

1.07107 

0.330331 

0.459301 

0.471418 

0.0125 

1.1427 

0.356006 

0.459289 

0.465441 

0.00625 

1.28276 

0.396078 

0.459294 

0.462488 

0.001 

9.14247 

0.635545 

0.459788 

0.459916 

Table  4.1:  The  intermediate  states  and  shock  speed  solutions  from  Equation 
(3.37)  based  on  the  precursor  thickness  b.  As  the  precursor  decreases,  both  hi 
and  fa  increase  and  the  shock  speeds  converge. 

The  intermediate  film  thickness  hl  and  particle  concentration  0,:  increase  as  the 
height  of  the  precursor  b  decreases.  For  the  shock  speeds,  a  smaller  precursor  leads 
to  the  trailing  shock  speed  Si  staying  relatively  the  same,  but  the  leading  shock 
speed  S2  slows  down  and  approaches  si .  These  results  agree  with  the  previous 
ones  related  to  solving  the  system  of  conservation  laws  [CBH08,  ZDB05].  For 
this  model,  the  smallest  precursor  for  which  a  solution  exists  is  b  «  9  x  10-4 
[CBH08].  A  precursor  close  to  this  case,  b  =  0.001,  produces  shocks  speeds  which 
are  close  together  and  an  intermediate  particle  concentration  near  the  maximum 
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packing  fraction.  An  alternative  settling  function  that  permits  solutions  with 
smaller  precursors,  /b(0)  =  (1  —  0 / <j>max)5 ,  is  examined  in  Cook  et  al.  [CBH08]. 

The  boundary  conditions  for  h  are  Dirichlet  in  front  and  behind,  in  the  re¬ 
direction,  the  flow  and  Neumann  on  the  sides,  in  the  ^/-direction.  The  same  is 
employed  for  <f>.  In  addition,  all  third  derivatives  in  h,  normal  to  the  boundary, 
are  set  to  0.  More  specifically,  for  a  rectangular  domain  with  length  X0  and  width 
Y0,  the  boundary  conditions  are 


^(0;2/)  I;  hxxa :(0)2/)  0)  h(Xo ,  r/)  bf  hxxx(Xo,  y')  0, 
h y  (x,  0)  0,  foyyy  (+  0)  0?  Ay  (x ,  4  o)  d  ,  ^yyy  ^  o)  d, 

0(0 ,y)  =  00,  0(^0 ,y)  =  00,  <t>y{x,  0)  =  0,  ( Py(x,Y0 )  =  0. 

The  simulations  are  all  run  using  moving  reference  frames,  with  the  speed  of  the 
frame  determined  as  in  Section  3.4. 

The  code  is  written  in  parallel  using  the  C++  OpenMP  package.  This  choice 
of  parallelization  was  made  since  the  majority  of  calculations  are  done  via  for 
loops  and  OpenMP  works  well  with  loop-heavy  code.  This  includes  the  calcula¬ 
tion  of  all  finite  differences  and  the  solves  along  rows  and  columns  associated  with 
the  ADI  part  of  the  scheme.  This  is  especially  useful  since  rows/columns  can  be 
solved  independently  of  each  other  for  each  equation.  In  addition,  writing  special 
solvers  for  linear  systems  of  equations  across  multiple  processors  [NS09,  PS06]  is 
avoided  by  this  approach.  The  speed-up  attained  using  N  processors  is  calcu¬ 
lated  by  dividing  the  runtime  for  one  processor  by  the  runtime  for  N  processors 
(Speed-Up  =  Time(l  Processor)/Time(iV  Processors)). 

Based  on  Figure  4.1,  the  scaling  is  close  to  linear  up  to  4  processors,  with 
a  small  drop-off  in  performance  as  the  number  increases.  This  almost-linear 
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Figure  4.1:  The  speed-up  gained  by  going  from  1  to  IV  processors  using  OpenMP. 
The  line  y  =  IV  is  shown  as  a  point  of  reference. 

behavior  is  a  result  of  all  of  the  code,  outside  of  a  few  minor  calculations  and  the 
recording  of  the  data,  being  amenable  to  parallelization. 

To  test  some  preferences  that  need  be  chosen  a  priori  in  the  simulation,  we 
conducted  short-time  tests  to  gauge  the  effectiveness  of  each  approach.  The  ones 
considered  here  are  (a)  whether  to  time-lag  or  extrapolate  the  approximate  terms 
and  (b)  whether  or  not  to  perform  iterations  past  a  single  solve  to  improve  the 
approximate  terms,  and  therefore  the  solution  at  each  timestep  (see  Table  4.2). 


(a)  Approximate  Terms 

Time-Lagged 

Extrapolation 

(b)  Number  of  Iterations 

One  Iteration 

Iterations 

Table  4.2:  The  two  choices  to  be  made  when  implementing  the  numerical  scheme. 
One  must  choose  whether  to  (a)  time-lag  or  extrapolate  the  approximate  terms 
and  (b)  whether  or  not  to  perform  additional  iterations  past  the  initial  solve. 
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Consider  an  initial  condition  of  </>0  =  0.3  and  a  front  perturbed  from  Riemann 
initial  data,  h(x,y,  0)  =  1  far  behind  the  front,  h(x,y,  0)  =  0.05  far  ahead  of  the 
front.  At  the  jump  from  fluid  to  precursor,  the  shape  of  the  front  is  given  as 
Xfront  =  X0/2  —  cos(2ny / Ya) .  This  initial  data  is  then  smoothed  using  hyperbolic 
tangent  and  matched  to  the  boundary  condition  (see  Figure  4.4).  This  has  the 
effect  that  the  initial  timestep  can  be  taken  more  leniently. 

We  ran  this  initial  simulation  for  each  of  the  four  combinations  in  Table  4.2  to 
time  t  —  1  and  the  maximum  timestep  allowed,  average  number  of  iterations  per 
timestep,  and  the  total  runtime,  in  seconds,  are  listed  in  the  table  below  (Table 
4.3).  This  and  Table  4.4  provide  some  global  measures  to  compare  the  different 
schemes  rather  than  illustrating  convergence  studies  for  any  particular  method. 
The  choice  of  t  =  1  was  made  as  the  timestep  changes  dramatically  over  this  time 
interval  and  can  provide  insight  as  to  which  methods  seem  practical  for  long-time 
runs.  Since  adaptive  timestepping  is  utilized  here,  the  tolerances  are  tuned  so  as 
to  ensure  that  the  simulation  stays  stable,  not  only  to  time  t  —  1  but  for  some 
time  afterwards  as  well  (it  is  taken  up  to  t  —  100  in  this  case,  which  is  the  length 
of  the  long-run  simulations). 


^^max 

Avg.  Iter. 

Runtime 

Time-Lagged  and  One  Iteration 

0.000568341 

1.0 

518.2 

Time-Lagged  and  Iterations 

0.00183296 

2.20997 

601.468 

Extrapolation  and  One  Iteration 

4.07743  x  10~5 

1.0 

19596.1 

Extrapolation  and  Iterations 

0.00486338 

1.29668 

376.603 

Table  4.3:  Results  for  time  t  —  1  based  on  various  choices  for  implementation. 


Using  Iterations  performs  well  for  both  choices  of  approximate  terms  in  that 
the  total  runtimes  are  low,  the  maximum  timesteps  are  large,  and  the  number 
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of  iterations  stays  close  to  1.  Between  these  two,  Extrapolation  and  Iterations 
does  best,  with  nearly  one  fewer  iteration  required  per  timestep,  on  average, 
and  a  runtime  that  is  37%  shorter.  Performing  One  Iteration ,  the  runtime  for 
Time-Lagged  is  in  between  the  two  cases  with  Iterations ,  but  for  Extrapolation , 
it  performs  poorly,  producing  a  runtime  that  is  33  to  52  times  worse  than  the 
other  three  options.  This  is  due  to  the  small  maximum  timestep  that  is  associated 
with  this  approach,  which  is  14  to  119  times  smaller  than  the  other  three.  At  this 
point,  it  makes  sense  to  discard  the  Extrapolation  and  One  Iteration  approach 
due  to  its  excessive  runtime  and  explore  the  remaining  ones. 

Under  the  same  conditions,  we  ran  a  longer  simulation,  this  time  to  t  =  100. 
Using  the  best  remaining  options,  we  can  glean  some  idea  as  to  which  one(s)  will 
work  best  for  a  longer  simulation. 


^^max 

Avg.  Iter. 

Runtime 

Time-Lagged  and  One  Iteration 

0.00107169 

1.0 

17811.3 

Time-Lagged  and  Iterations 

0.00329173 

2.95498 

13153.8 

Extrapolation  and  Iterations 

0.0106161 

2.01204 

3364.93 

Table  4.4:  Results  for  time  t  =  100. 


Comparing  Tables  4.3  and  4.4,  the  maximum  timestep  for  each  approach 
has  increased.  Using  Iterations ,  the  average  number  has  gone  up  for  both  Time- 
Lagged  and  Extrapolation.  However,  the  average  number  of  iterations  per  timestep 
for  Extrapolation  is  approximately  one  iteration  fewer  than  for  Time- Lagged.  Also 
the  runtime  takes  about  2.9  times  longer  for  Time-Lagged  compared  to  Extrapo¬ 
lation.  One  can  see  the  benefit  of  performing  iterations  instead  of  using  a  smaller 
timestep  in  comparing  the  results  for  Time-Lagged  and  One  Iteration  and  Time- 
Lagged  and  Iterations.  Time-Lagged  and  One  Iteration  advances  the  solution 
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approximately  the  same  time  forward  with  three  timesteps  as  Time-Lagged  and 
Iterations  does  with  one  timestep  and  three  iterations.  However,  doing  two  extra 
timesteps  costs  more  than  two  extra  iterations,  as  seen  in  their  respective  run¬ 
times.  This  is  because  the  explicit  terms  do  not  need  to  be  re-calculated  for  each 
iteration  while  they  do  for  each  timestep.  Therefore,  the  only  two  options  which 
make  sense  to  use  are  the  ones  involving  Iterations.  Of  these,  Extrapolation  is 
the  clear  favorite. 
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Figure  4.2:  The  adaptive  timestep  up  to  time  t  =  20.  The  timestep,  At,  is 
recorded  in  intervals  of  0.25  for  the  three  cases.  Extrapolation  and  Iterations 
has  a  significantly  larger  timestep  than  either  Time-Lagged  and  One  Iteration  or 
Time-Lagged  and  Iterations. 

In  Figure  4.2,  we  see  that  by  time  t  —  8,  all  three  approaches  have  settled  into 
a  respective  timestep.  The  timestep  for  Extrapolation  and  Iterations  does  best, 
followed  by  Time-Lagged  and  Iterations  and  Time-Lagged  and  One  Iteration.  The 
timestep  for  Extrapolation  and  Iterations  is  3.2  times  better  than  Time-Lagged 
and  Iterations  and  9.9  times  better  than  Time-Lagged  and  One  Iteration.  The 


-Time-Lagged  and  One  Iteration 
Time-Lagged  and  Iterations 
Extrapolation  and  Iterations 
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benefit  of  the  larger  timestep  for  both  approaches  with  Iterations  is  partially 
offset  by  the  need  for  extra  calculations  related  to  the  iterations. 


0  20  40  60  80  100 

Time  (t) 

Figure  4.3:  The  number  of  iterations  up  to  time  t  =  100.  The  iterations  are 
recorded  in  intervals  of  0.25  for  the  two  cases.  Using  Extrapolation  and  Itera¬ 
tions  does  better  than  Time-Lagged  and  Iterations  in  terms  of  fewest  number  of 
iterations. 

Figure  4.3  shows  the  number  of  iterations  required  throughout  the  simulation. 
For  Extrapolation  and  Iterations ,  the  increase  in  iterations  approximately  between 
times  t  —  20  and  t  =  30  corresponds  to  the  finger  forming  and  stretching  out 
ahead  of  the  flow  in  the  film  thickness  and  the  particle-rich  ridge  growing  higher 
and  outlining  the  finger.  While  the  number  of  iterations  jumps  once  to  3  and 
then  back  down  to  2  for  Extrapolation  and  Iterations ,  it  remains  constant  at  3 
for  Time-Lagged  and  Iterations.  The  cost  of  storing  extra  data  and  performing  a 
small  computation  to  find  the  extrapolated  approximations  seems  a  small  price 
to  pay  to  save  one  iteration  per  timestep,  which  includes  recalculating  values 
involving  the  approximate  terms  and  performing  the  ADI  solves. 

Using  the  simulation  data  up  to  t  —  100,  we  can  examine  the  effects  of  the 
initial  perturbation  graphically.  For  the  film  thickness,  a  small  capillary  ridge 
forms  in  the  center  of  the  perturbation  (Figure  4.5)  and  begins  to  stretch  out 
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Figure  4.4:  The  initial  film  thickness.  It  is  perturbed  by  a  cosine  wave  along  y 
and  smoothed  along  x  by  hyperbolic  tangent. 


Figure  4.5:  Film  thickness  (left)  and  particle  concentration  (right)  at  time  t  =  25. 
A  small  ridge  forms  in  both,  with  the  highest  point  in  the  perturbation. 

ahead  of  the  bulk  flow  (Figures  4.6  and  4.7).  This  is  the  well-known  fingering 
instability  present  in  thin  film  flows.  For  the  particle  concentration,  a  particle- 
rich  ridge  initially  forms  at  the  contact  line  (Figure  4.5)  and,  as  the  fingering 
instability  evolves,  outlines  the  shape  of  the  finger  (Figures  4.6  and  4.7).  Directly 
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Figure  4.6:  Film  thickness  (left)  and  particle  concentration  (right)  at  time  t  =  50. 
A  fingering  instability  and  particle- rich  ridge  form. 


Figure  4.7:  Film  thickness  (left)  and  particle  concentration  (right)  at  time 
t  =  100.  The  fluid  finger  stretches  out  ahead  of  the  bulk  flow.  The  particle-rich 
ridge  increases  in  concentration  and  has  a  higher  concentration  in  and  around 
the  fingering  instability. 

behind  the  ridge,  a  pocket  of  lower  concentration  forms.  The  interior  of  the 
finger  is  slowly  encroached  upon  by  the  particles  that  have  accumulated  near  the 
back  and  sides  of  the  finger.  This  can  be  seen  in  Figure  4.7  as  an  interior  layer 
along  the  inside  of  the  particle-rich  ridge.  It  is  possible  that  this  phenomenon 
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is  not  physical,  meaning  that  it  occurs  only  in  the  simulations  and  not  in  the 
experiments,  and  may  be  a  result  of  the  current  model  not  containing  all  of  the 
necessary  physics. 

4.2  Comparison  to  Experiments 

Experiments  for  constant-volume  particle-laden  thin  film  flows  have  been  com¬ 
pared  in  one  dimension  to  the  solution,  both  analytically  and  numerically,  for 
constant-volume  clear  fluid  flows.  The  average  front  position  for  clear  fluids  is 
given  by  a  power  law,  where  the  location  of  the  front  scales  like  Ct 1//3,  where 
C  is  a  scaling  constant  [Hup82],  Ward  et  al.  [WWG09]  compare  the  average 
front  position  of  the  flow  to  this  scaling  and  find  agreement  for  particle  concen¬ 
trations  below  (p  —  0.45  and  deviations  at  later  times  for  higher  concentrations. 
Grunewald  et  al.  [GLM10]  compare  the  average  front  position  to  a  re-derived 
one-dimensional  model,  based  on  results  from  Huppert  [Hup82]  with  a  precur¬ 
sor,  and  to  experiments  and  numerical  solutions  of  the  one-dimensional  problem. 
The  Ct1/3  scaling  appears  valid  for  concentrations  of  0.25  to  0.45,  and  the  scaling 
constant  for  experiments  and  numerics  are  within  20%  of  the  theoretical  constant 
(see  Chapter  5  for  a  full  discussion).  We  seek  to  compare  the  numerical  solution 
in  two  dimensions  to  images  of  experiments,  taking  into  account  that  variations 
occur  across  the  front  of  the  flow. 

We  use  1000  cSt  polydimethylsiloxane  (PDMS),  a  silicone  oil,  for  the  liquid 
component  of  the  fluid.  For  the  particles,  glass  beads  with  diameters  in  the  range 
of  250  —  425  /jrn  are  used.  The  two  components  are  then  well-mixed  and  released 
down  an  inclined  plane  from  a  reservoir.  This  corresponds  to  a  constant-volume 
experiment,  whereas  our  simulations  are  constant-flux.  The  approximation  of  a 
constant-volume  problem  by  a  constant-flux  one  may  be  invalid  at  early  times, 
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as  the  height  of  the  fluid  will  be  changing  rapidly.  However,  the  height  of  the 
flow  changes  slower  at  later  times,  at  which  point  a  constant-flux  approximation 
may  be  valid. 

The  experiment,  which  we  will  compare  to  simulation,  is  a  fluid  of  approxi¬ 
mately  90  cm3  containing  a  volume  which  is  35.9%  particles.  The  plane  is  inclined 
at  a  32-degree  angle.  The  fluid  is  allowed  to  flow  down  the  plane,  which  is  14  cm 
across  and  90  cm  down.  In  the  experiments,  the  flow  starts  out  close  to  uniform 
across  the  front,  away  from  the  edges,  and  over  time  develops  instabilities,  in 
the  form  of  fingers  stretching  out  ahead  of  the  bulk  flow.  Since,  for  simulations, 
starting  with  a  uniform  front  along  the  ^/-direction  leads  to  a  uniform  solution, 
we  start  the  simulation  some  time  after  the  start-time  to  add  a  perturbation  to 
the  initial  data,  which  induces  the  type  of  behavior  seen  in  the  latter  stages  of 
the  experiments. 

In  order  to  avoid  simulating  the  problem  over  the  entire  domain,  we  truncate 
the  solution  near  the  front  and  treat  the  problem  locally  as  being  constant-flux. 
We  are  interested  in  the  dynamics  of  finger  formation  during  which  time  the  film 
thickness  only  changes  by  at  most  20%  ,  so  a  local  approximation  by  constant-flux 
is  reasonable. 

We  use  two  images,  taken  three  minutes  apart,  to  compare  with  the  simulation 
(Figures  4.8  and  4.10).  The  first  image  is  taken  when  the  front  of  the  flow  has 
reached  approximately  53  cm  down  the  plane.  The  shape  of  the  front  is  parabolic- 
like  with  two  large  perturbations  at  either  end  of  the  front.  In  between,  smaller 
perturbations  exist  which  lead  to  fingering  instabilities.  The  two  outer  perturba¬ 
tions  lead  to  longer  and  thicker  fingers  than  the  smaller  inner  perturbations.  We 
take  a  front  similar  to  this  in  our  simulation. 

The  scales  for  a  constant-flux  problem  can  be  taken  from  Cook  et  al.  [CBH08], 
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Figure  4.8:  The  initial  condition  of  the  experiment,  used  for  comparing  with  the 
simulation.  At  this  point,  the  front  of  the  flow  has  begun  to  develop  perturba¬ 
tions,  which  will  lead  to  fingering  instabilities. 

which  are  the  same  as  for  the  clear  fluid  case.  The  height  scale  is  taken  to  be 
/to  —  1  mm.  The  length  scale  is  Xo  =  (Z2/io) 1//3 ,  where  the  capillary  length,  /, 
is  l  =  The  constants  are  7,  the  coefficient  of  surface  tension;  pi,  the 

liquid  density;  and  g\\,  the  component  of  gravity  parallel  to  the  inclined  plane. 
The  time  scale  is  t0  =  (3pi/^f)xol2 /h^,  where  /q  is  the  dynamic  liquid  viscosity. 
The  capillary  number  is  given  by  Ca  =  pix0/^/t0  =  /tg/3/2. 

The  scales,  given  these  parameters,  are  /to  =  0.001  m,  xq  =  0.00161396  m, 
/  =  0.00205041  m,  t0  =  0.93235  s,  and  Ca  =  0.0792863.  Using  this,  we  can 
construct  an  initial  condition  which  resembles  the  experiment  and  will  produce 
similar  results.  This  is  done  by  measuring  the  features  of  the  initial  image  and 
creating  a  similar  condition.  While  the  flow  in  the  experiment  is  asymmetric,  we 
take  a  symmetric  initial  condition  in  the  simulation  which  has  features  that  are 
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approximately,  in  both  location  and  size,  the  same  as  in  the  experiment.  The 
track  is  taken  to  be  86.75  units  wide  (rounded  up  to  the  nearest  0.05  increment, 
which  is  the  value  of  Ax,  Ay),  which  corresponds  to  the  14  cm  wide  track.  The 
precursor  in  the  simulation  is  set  to  b  =  0.05,  as  in  the  previous  simulations. 

A  moving  reference  frame  is  used  since  this  is  taken  to  be  a  constant-flux 
problem  locally.  The  speed  of  the  moving  reference  frame  is  approximately 
s  =  0.343198,  calculated  as  in  Section  3.4.  Running  a  simulation  over  the  course 
of  three  minutes  leads  to  a  distance  traveled  for  the  frame  of  approximately 
10.69  cm,  where  the  actual  displacement,  based  on  the  experiment,  is  around 
12  cm,  so  using  the  constant-flux  assumption  seems  to  produce  a  decent  approx¬ 
imation  of  the  distance  the  fluid  will  flow. 
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Figure  4.9:  The  initial  condition  for  the  film  thickness,  h,  used  in  the  simulation. 
This  is  an  artificially-created  starting  condition  to  be  representative  of  the  state 
shown  for  the  experiment.  The  height  is  in  mm. 

The  initial  data  is  generated  using  a  sine  wave  to  form  the  two  large  perturba- 
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tions  and  the  space  away  from  the  edges.  The  three  fingers  that  develop  between 
these  two  perturbations  are  simulated  with  a  cosine  wave  of  small  amplitude, 
0.25  in  dimensionless  units  (Figure  4.9).  The  simulation  is  run  to  t  —  193.06,  the 
equivalent  of  three  minutes  of  real-time. 


Figure  4.10:  The  evolution  of  the  experiment  after  three  minutes.  The  fingering 
instability  starts  to  form  at  the  front. 

Over  the  course  of  the  three  minutes,  the  exterior  of  the  outer  fingers  in  the 
experiment  go  from  4  cm  and  6.5  cm  on  the  left  and  right,  respectively,  to  7.5  cm 
and  12  cm.  The  interior  of  these  fingers  go  from  less  than  1  cm  on  each  side  to 
about  3  cm.  The  interior  fingers  are  not  discernable  in  the  initial  image.  The  flow 
as  a  whole,  measured  from  where  the  fluid  touches  the  walls,  has  moved  about 
11  cm  down  the  plane.  The  interior  fingers  in  the  experiment  extend  about  0.5  cm 
ahead  of  the  flow. 

In  the  simulation  (Figure  4.11),  the  moving  reference  frame  accounts  for 
10.69  cm  of  movement,  so  the  position  where  the  fluid  touches  the  edges  has 
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Figure  4.11:  The  evolution  of  the  him  thickness,  h,  in  the  simulation  after  three 
minutes.  Both  the  experiment  and  simulation  exhibit  a  fingering  instability,  but 
the  instability  in  the  simulation  is  less  pronounced.  The  height  is  in  mm. 
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Figure  4.12:  Particle  concentration,  </>,  for  the  him  thickness  in  Figure  4.11. 
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moved  approximately  11.4  cm.  The  evolution  of  the  fingers  in  the  simulation  is 
slightly  less  pronounced  than  in  the  experiment.  This  is  likely  dne  to  the  sim¬ 
ulation  initially  undergoing  a  transient  state  where  the  fluid  travels  slower  than 
at  later  times,  while  the  transient  in  the  experiment  has  occurred  prior  to  this 
three-minute  interval.  The  exterior  of  the  outer  fingers  is  approximately  4.2  cm 
and  interior  1.2  cm.  The  interior  fingers  extend  ahead  of  the  flow  about  0.8  cm. 
The  tip  of  the  longest  finger  in  the  experiment  has  moved  15  cm  while  in  the 
simulations,  it  has  advanced  approximately  11.4  cm.  The  tips  of  the  fingers,  in 
the  ^-direction,  reach  up  to  1.37  mm. 

The  particle  concentration  cannot  be  determined  accurately  at  the  particle- 
rich  ridge  in  the  experiment,  but  the  increased  opacity  at  the  leading  edge  of  the 
flow  indicates  an  increase  in  the  concentration,  relative  to  the  ambient  concen¬ 
tration.  This  change  in  shade  is  approximately  1  cm  long  in  the  direction  of  the 
flow.  In  the  simulation  (Figure  4.12),  the  thickness  of  the  ridge  ranges  from  0.6 
to  1.1  cm,  which  is  consistent  with  the  experiment. 

4.3  Summary 

Schemes  originally  derived  for  numerically  solving  high-order  parabolic  problems 
have  recently  been  extended  to  high-order  systems,  such  as  the  case  of  surfactants 
and  particle-laden  thin  films.  Handling  the  higher-order  terms  in  a  practical  way 
is  necessary  for  fast  and  efficient  computation.  The  scheme  we  have  presented 
in  Chapter  3  for  particle-laden  thin  film  flow  provides  an  easy-to-program  and 
effective  way  to  solve  this  high-order  coupled  system.  This  scheme  can  provide  a 
blueprint  for  approaches  to  solving  similar  problems. 

The  numerical  scheme  developed  for  particle-laden  thin  film  flow  has  several 
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nice  attributes.  The  timestep  required  for  this  scheme  is  in  the  range  of  0( Ax2), 
which  is  much  better  than  the  0( Ax4)  for  a  fully  explicit  scheme.  The  structure 
of  the  scheme  allows  for  the  possibility  of  solving  each  equation  with  its  own 
unique  timestep  for  better  efficiency,  as  the  particle  concentration  is  typically  the 
equation  that  fails  the  timestep  restriction  criteria.  The  linear  algebra  problem 
that  results  from  an  implicit  time  discretization  along  with  the  nonlinearity  is 
reduced  to  a  series  of  tri-  and  pentadiagonal  solves,  which  can  be  done  in  parallel 
along  the  rows/columns  of  the  grid. 

The  parallelization  of  the  code  is  straightforward  using  OpenMP.  The  loops 
for  computing  the  explicit  and  approximate  terms  as  well  as  the  solves  along 
rows  and  columns  can  be  done  in  parallel,  leading  to  a  code  that  scales  close  to 
linearly  for  up  to  8  processors,  getting  close  to  8  times  speed-up.  Adding  OpenMP 
implementation  to  C++  code  on  any  multicore  machine  is  easy  to  implement, 
as  it  only  requires  adding  a  few  lines  of  code  to  existing  for  loops  and  needs  no 
managing  of  the  movement  of  data  on  the  programmer’s  part.  Since  the  code  is 
predominantly  such  loops,  it  is  easy  to  parallelize  and  is  highly  effective  in  getting 
better  runtimes. 

Implementing  Iterations  within  each  timestep,  which  is  first  presented  in  Wi- 
telski  and  Bowen  [WB03],  but  not  used  in  Warner  et  al.  [WCM05],  seems  to  work 
best  for  this  problem,  in  terms  of  allowing  for  a  larger  timestep  and  producing  an 
accurate  solution.  Among  the  choices  for  the  approximate  terms  when  performing 
Iterations ,  Extrapolation  seems  to  produce  the  best  runtime  and  fewest  iterations. 
Implementation  requires  only  storing  an  extra  set  of  data  used  in  extrapolating 
the  approximate  terms  but,  using  the  adaptive  timestepping  discussed  here,  this 
data  is  stored  anyway. 

The  choice  of  Extrapolation  and  Iterations  may  work  best  for  this  problem, 
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but  for  other  problems  or  initial  conditions,  another  choice  may  fare  better.  It 
is  recommended,  as  in  this  case,  that  a  short-term  simulation  be  performed  for 
the  different  choices  of  approximate  terms  and  whether  or  not  to  perform  extra 
iterations.  The  small  cost  of  these  short  runs  may  allow  for  a  more  efficient  run 
for  actual  simulations.  It  is  also  recommended  that  one  examines  the  results 
to  make  sure  that  the  scheme  is  not  only  fast  with  the  choice,  but  sufficiently 
accurate. 

The  numerical  solution  agrees  reasonably  well  with  the  behavior  seen  in  ex¬ 
periments.  This  is  in  part  because  the  model  was  derived  for  the  case  when  a 
particle-rich  ridge  forms.  This  is  seen  in  the  experiments  for  high  angles  of  incli¬ 
nation  and  high  concentrations,  but  will  occur  in  the  model  for  all  concentrations 
and  angles.  The  particle-rich  ridge  in  the  simulations  is  two  thin  layers  of  parti¬ 
cles,  one  which  originates  at  the  front  of  the  flow  and  the  other  from  the  troughs 
of  the  emerging  fingers,  which  may  not  be  physical. 

The  current  model  assumes  a  constant,  or  average,  particle  concentration 
throughout  the  fluid  layer  in  the  ^-direction.  The  same  is  true  for  the  velocity, 
which  is  averaged  in  the  ^-direction.  Theory  exists  for  the  vertical  movement 
of  the  particles  [C00O8],  whether  they  will  settle  to  the  inclined  plane  or  form  a 
ridge,  and  incorporating  this  behavior  into  a  new  model  is  the  current  research 
of  the  authors.  It  is  hoped  that  the  current  numerical  scheme  will  be  adaptable 
to  this  new  model. 
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CHAPTER  5 


Self-Similarity  of  Particle-Laden  Flow  at 
Constant  Volume 

This  chapter  compares  the  behavior  of  a  constant-volume  particle- laden  fluid  flow 
down  an  incline  with  a  clear  fluid  flow.  We  compare  the  analytical  solutions  of 
the  clear  fluid  flow  to  both  physical  experiments  and  numerical  simulations  of 
the  particle-laden  fluid  flow,  using  a  model  proposed  in  Cook  et  ah  [CBH08]  and 
Zhou  et  al.  [ZDB05]. 

In  Huppert  [Hup82] ,  a  simple  scaling  law  is  derived  for  the  average  front  posi¬ 
tion  x(t)  ~  Cf1//3  in  the  case  of  clear  fluids.  Comparison  to  both  the  particle-laden 
flow  lubrication  model  and  to  physical  experiments  suggests  that  the  Huppert 
scaling  law  is  still  valid  to  leading  order  for  particle-laden  fluids  with  moderate 
particle  concentrations  in  the  range  25-45%.  In  this  range,  the  particle-laden  fluid 
still  behaves  fluid-like,  and  settling  of  the  particles  is  present  but  does  not  domi¬ 
nate  the  large-scale  dynamics.  The  effects  of  settling  in  the  direction  of  the  flow 
can  be  visually  observed  in  the  experiments  as  a  particle-rich  ridge  at  the  leading 
edge  of  the  flow.  We  compare  different  settling  functions  in  our  model  to  analyze 
this  effect  numerically.  We  also  note  that  the  lubrication  models  with  settling 
require  a  precursor;  they  are  singular  at  vanishing  precursor  [CBH08,  ZDB05]. 
Thus,  it  makes  sense  to  compare  the  dynamics  of  the  lubrication  model  with 
settling  to  an  exact  solution  of  the  problem  without  settling  and  with  precursor. 
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5.1  One-Dimensional  Lubrication  Model 


The  flow  of  a  particle-laden  fluid  down  an  inclined  plane  has  been  modeled  by  a 
system  of  conservation  laws  [CBH08,  Coo07], 


(5.1) 

{</>h)t  +  {(/)hvp)x  =  0. 

(5.2) 

This  is  similar  to  the  system  of  conservation  laws  (3.33)-(3.34)  in  Section  3.4. 
The  total  volume  of  the  fluid  is  conserved  in  (5.1),  where  h(x,t)  is  the  height 
of  the  liquid-particle  mixture  at  position  x  (oriented  down  the  inclined  plane) 
and  time  t.  We  assume  uniformity  in  the  transverse  direction.  The  total  volume 
of  particles  <ph  is  conserved  in  (5.2),  where  0  is  the  particle  concentration  in  the 
fluid.  The  average  velocity  of  the  fluid  utot  is  a  depth-average  as  well  as  a  volume- 
average  of  the  speeds  of  the  liquid  and  the  particles.  The  velocity  of  the  particles 
Up  consists  of  utot  and  an  extra  term  (1  —  0)urei  due  to  the  sedimentation  of  the 
more  dense  particles  in  the  fluid: 

vp  =  vtot  +  (1  -  0)ure i. 

We  assume  that  the  particle  concentration  0  is  uniform  in  the  direction  perpen¬ 
dicular  to  the  inclined  plane.  It  is  shown  in  Cook  [C00O8]  that  this  assumption 
can  be  improved  to  a  distribution  that  is  stationary  in  time  which  balances  sedi¬ 
mentation  and  shear-induced  migration.  We  refer  to  Cook  [C00O8]  for  a  precise 
discussion.  The  depth-average  of  the  velocity  of  the  fluid,  utot  found  in  (5.1),  can 
then  be  derived  via  standard  techniques  in  lubrication  theory  for  thin  liquid  films 
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[Coo07],  which  is  valid  for  fluids  with  small  Reynolds  numbers  and  characteristic 
height  much  smaller  than  the  characteristic  length.  The  first-order  term  describes 
its  dominant  behavior  as: 


<5'3> 

Note  that  the  terms  discussed  in  this  chapter  (e.g.,  p(^) ,  p((f>) , . . .)  have  not  been 
nondimensionalized.  Following  the  modeling  suggested  in  Krieger  [Kri72]  and 
Stickel  and  Powell  [SP05],  we  model  the  viscosity  with  the  empirically  derived 
model 


P-{4>)  =  ni(  1  -  0/0max)  2  (5.4) 

with  dynamic  liquid  viscosity  pi.  The  value  0max  is  as  defined  in  Chapter  2.  We 
use  0max  =  0.58,  which  is  our  experimentally  determined  value  [WWG09].  This 
value  is  also  used  in  the  literature;  see  e.g.,  [ZHL00].  The  density  is  a  linear 
combination  of  the  density  of  the  liquid  pi  and  the  density  of  the  particles  pp 

p{<t>)  =Pi(  1  -  4>)  +  PP(p- 

The  other  constant  in  (5.3)  is  the  component  of  gravitational  acceleration  g\\  = 
\g\  since,  where  a  is  the  inclination  angle  of  the  plane. 

The  velocity  up  in  the  second  conservation  law  (5.2)  requires  a  more  explicit 
description,  since  the  theory  for  particle-laden  flow  is  still  relatively  novel  and  con¬ 
tinues  to  present  many  open  questions,  especially  for  shear-driven  flows.  Recall 
that  the  lubrication  approximation  used  to  derive  (5.3)  employs  a  depth-averaged 
velocity.  We  also  employ  a  depth-averaged  model  for  urei  as  in  Chapter  2,  which 
we  assume  is  a  product  of  three  factors: 


53 


vtei  =  VJ(<f})w(h). 


(5.5) 


The  Stokes  settling  velocity  Vs  is  given  as: 

T f  _  2  (Pp  ~  Pl)g\\a2 

*s  —  n 

9  Hi 

The  other  factors  account  for  phenomena  that  reduce  the  speed  of  a  single  par¬ 
ticle:  hindered  settling  from  adjacent  particles  and  slowing  due  to  proximity  of 
particles  to  the  substrate. 

A  classical  model  for  hindered  settling  was  proposed  by  Richardson  and  Zaki 
[RZ54a]  and  Buscall  et  al.  [BG082]: 

f{4>)  =  (1  -  </>/0max)m  (5.6) 

with  0max  =  1  and  empirically  determined  exponent  m  =  5.1.  Cook  [Coo07]  mod¬ 
ified  the  function  to  include  the  maximum  packing  fraction  of  particles  (0max  = 
0.67  in  their  case)  to  avoid  singular  shocks  in  solutions  to  the  Ricmann  problem 
for  (5.1)  and  (5.2)  that  occur  when  0max  =  1-  This  form  also  ensures  that  sed¬ 
imentation  stops  once  the  maximum  concentration  is  reached.  We  will  compare 
results  for  m  —  1  and  m  —  5  to  probe  the  effect  of  the  exponent  in  the  hin¬ 
dered  settling  function.  This  is  particularly  relevant  when  comparing  numerical 
results  and  physical  experiments,  since  the  division  by  the  maximum  packing 
fraction  may  have  altered  the  appropriate  choice  of  exponent  for  comparison  to 
experiments.  Note  that,  although  the  singular  limits  for  both  functions 

lim  f (</>)  =  0,  lim/(0)  =  l, 

<P  ^0max  0—^0 
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are  appropriate,  we  will  not  consider  extreme  values  of  0,  since  the  comparison 
to  experiments  is  most  appropriate  for  moderate  concentrations  of  particles. 

The  third  factor  in  (5.5)  models  the  slowing  of  particles  due  to  their  proximity 
to  the  substrate,  called  the  wall  effects  function  [ZDB05,  Coo07]  (see  Chapter 
2).  The  full  system  of  equations  for  h(x,t )  and  4>(x,t)  is  now  fully  specified  by 
incorporating  (5.3)  into  (5.1)  and  (5.5)  into  (5.2). 

As  has  been  shown  for  clear  fluids  [BB97],  a  first-order  model  such  as  the 
one  proposed  here  can  correctly  capture  quantities  such  as  front  speed  but  does 
not  contain  the  physics  necessary  to  model  fingering  due  to  surface  tension  and 
the  component  of  gravity  normal  to  the  substrate.  Numerical  evidence  in  Cook 
[Coo07]  for  constant-flux  boundary  conditions  indicates  that  higher-order  terms 
smooth  solutions  but  do  not  affect  the  speed  of  the  leading  front  of  the  him. 
Nevertheless,  quite  a  lot  of  information  can  be  gained  by  studying  the  dominant 
physics  in  the  how  direction,  in  particular  with  regard  to  the  competition  between 
settling  of  particles  and  overall  motion  of  the  fluid. 

Nondimensionalizing  the  reduced  system  with  length  scale  ho  (half  of  the 
upstream  gate  height)  and  time  scale  t0  =  (sin  ct)f,  we  obtain 


ht  + 


( 9P{<t>)  ,3\ 

V  3/^(0)  ) 


0, 


(5.7) 


-  0)/(0)M/0j  =  0,  (5.8) 

which  we  solve  numerically  and  compare  to  experimental  results.  This  nondi- 
mensionalization  is  slightly  different  than  the  one  given  for  the  two-dimensional 
system  in  Chapter  2.  Note  that  this  system  has  been  studied  for  constant-hux 
boundary  conditions  [Coo07].  Since  the  physical  experiment  (described  in  Sec- 


(#)< 
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tion  5.3)  more  closely  resembles  a  constant  volume  of  particle-laden  fluid  than 
a  constant  flux  (for  the  purpose  of  numerical  simulations),  we  will  choose  initial 
conditions  to  approximate  a  constant  volume  of  fluid. 

5.2  Well-mixed  similarity  theory  including  precursor 

We  begin  our  discussion  of  similarity  solutions  by  considering  a  model  for  a  well- 
mixed  fluid  flow  down  an  inclined  plane.  Let  (po  be  the  initially  constant  particle 
concentration  in  the  fluid.  For  a  well-mixed  fluid  without  particle  settling,  the 
concentration  stays  constant  in  space  and  time.  We  compare  the  settling  model 
(5.7)  and  (5.8)  to  solutions  of  the  well-mixed  model: 

ht  +  ( hi)x  =  °>  (5-9) 

with 


p  =  p(0o)  and  p  =  p(<p0). 

This  model  from  Huppert  [Hup82]  captures  the  dominant  behavior  of  the  flow 
with  no  particles.  For  simulations  with  no  precursor,  the  solution  of  (5.9)  for  an 
initial  fluid  profile 


h(x,  0)  = 


/?,  0  <  x  <  L, 


0,  otherwise 

can  be  solved  explicitly.  We  consider  (5.9)  in  its  natural  time  scale 


T  =  f>t 

3 /Jo 
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After  an  initial  transient  time  r  >  r* 


3  L 

w 


a  similarity  solution 


h(x,  t ) 

develops  with 


\Jx/t,  0  <  x  <  xs(t), 
0,  otherwise 


x, 


(r)  =  Q/Ji)  '  t1'3  =  Ct1'3, 


(5.10) 


where 


c-i>re?r 


(5.11) 


The  constant  C,  when  measured  experimentally,  provides  a  method  for  measuring 
the  effective  viscosity  of  the  fluid  by  inverting  (5.11) 


b  =  | /32L2gp0C  3. 

The  effectiveness  of  this  quantity  as  a  measure  for  bulk  viscosity  of  the  particle- 
laden  fluid  is  explored  in  Ward  et  al.  [WWG09].  For  comparison  between  nu¬ 
merical  experiments  and  simulations,  we  focus  on  the  variation  in  the  scaling 
constant  C  for  experiments  with  and  without  particles. 

It  is  also  possible  to  consider  analytical  solutions  of  (5.9)  with  a  small  uniform 
precursor  layer  of  height  b.  We  consider  (5.9)  with  Riemann  initial  data 


h(x,  0) 


/3,  0  <  x  <  L, 

b,  otherwise  . 
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The  discontinuity  in  the  initial  data  at  x  —  0  will  immediately  become  a  rarefac¬ 
tion,  while  the  shock  at  x  —  L  will  persist.  At  a  critical  time 

L 

the  trailing  rarefaction  and  leading  shock  will  merge,  creating  a  wedge  shape  that 
continues  to  evolve  in  time  above  the  precursor.  The  shape  of  this  profile  can  be 
described  analytically  [Lax73,  Eva98].  Above  the  precursor  the  solution  still  has 
the  shape  of  a Jx/t.  The  excess  volume  above  the  precursor 


V  =  L(f3-b)  (5.12) 

remains  constant.  One  can  therefore  use  the  conservation  of  volume  to  determine 
the  speed  of  the  shock.  Let  xs  be  the  position  of  the  shock.  Then 


\J xjr  dx  —  b(xs  —  b2r)  =  V. 

Therefore  the  shock  position  is  defined  implicitly  by  the  solution  of 


<b2 1 


2 

3 


(xs\/xs/t  —  b3T^j  —b(xs  —  b2r)  —V  =  0.  (5.13) 

and  the  solution  of  the  double  Riemann  problem  with  precursor  after  rc  is 


( 


b, 


h{x,  t )  =  < 


X  <  b2r, 
b2r  <  X  <  xs, 


y  b,  x  >  xs. 


(5.14) 


To  explore  the  difference  in  evolution  of  the  solutions  analytically,  let 


H  =  \Jx/t 


(5.15) 


be  the  maximum  height  of  h.  Rewriting  (5.13)  in  V  and  H  gives: 

| H 3  -  V/t  =  bH 2  -  b3/ 3. 

This  illustrates  by  (5.15)  and  (5.12)  that,  as  long  as  b  <C  H,  the  right-hand 
side  is  negligible  and  xs  approximates  the  shock  position  of  the  solution  without 
precursor  (5.10).  Figure  5.1  shows  the  deviation  of  the  shock  positions.  The 
data  start  from  time  rc  after  which  (5.13)  is  valid.  For  increasing  time  and 
precursor,  the  relation  b  <C  H  gets  less  valid  and  therefore  the  deviation  of  the 
shock  positions  increases.  Note  that  for  the  settling  model  the  limit  b  — »  0  is 
singular  [ZDB05,  CBH08].  We  therefore  have  to  introduce  a  precursor  for  the 
numerical  simulations. 

5.3  Experimental  Results 

The  experimental  apparatus  consists  of  a  90  cm  long,  50  cm  wide  acrylic  sheet 
mounted  to  an  adjustable  stand  capable  of  inclination  angles  ranging  from  5°  to 
80°  (see  Figure  6.1).  Down  the  length  of  the  sheet  is  a  track  approximately  14 
cm  wide.  Side  walls,  approximately  3.2  cm  high  near  the  top  and  1.4  cm  near 
the  lower  part  of  the  track,  are  designed  so  the  fluid  does  not  escape  over  the 
sides.  Near  the  top  of  the  acrylic  sheet  is  a  gated  reservoir  from  which  a  finite 
well-mixed  volume  of  liquid  and  particles  is  released. 

The  experiments  shown  are  all  conducted  at  an  inclination  angle  of  45°  and 
particle  concentrations  of  25-45%,  in  increments  of  5%.  For  these  values,  we 
avoid  the  rapid  settling  of  the  particles  toward  the  substrate  associated  with  low 
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Figure  5.1:  Comparison  of  front  position  for  different  precursor  heights,  at 
L  =  2.5,  (3  =  2  and  0  =  0.4. 

concentrations  and  low  inclination  angles  which  leads  to  deposits  of  particles  at 
the  rear  of  the  flow  and  clear  fluid  fingers  at  the  leading  edge.  We  also  avoid 
particle  jamming,  clumping,  and  sliding  that  is  associated  with  higher  particle 
concentrations  and  higher  inclination  angles.  Additional  experimental  data  have 
been  collected  for  other  inclination  angles  [WWG09];  but  the  data  presented  here 
for  45°  provide  representative  results. 

The  particle-laden  fluid  is  a  mixture  of  1000  cSt  silicone  oil  (Clearco  Products) 
with  a  density  of  approximately  0.96  g/cm3,  and  soda-lime  glass  beads  (Ceroglass) 
with  a  density  of  approximately  2.5  g/cm3.  The  diameter  of  the  beads  is  0.025 
cm.  For  smaller  particles,  the  settling  in  the  direction  of  the  flow  to  form  a 
particle  ridge  at  the  front  is  less  evident. 

This  experimental  setup  has  a  Peclet  number  Pe  on  the  order  of  1010,  with 


60 


Pe  =  Vsa/D  where  Vs  =  (2/9 )(pp  —  pi)g\\a2 / pi  is  the  settling  velocity,  a  is  the 
particle  radius,  and  D  =  kiT/Qiipi  is  the  diffusivity  with  Boltzmann  constant 
kb  and  temperature  T.  We  note  that  a  high  Peclet  number  in  the  experiments 
suggests  that  we  may  neglect  Brownian  motion  in  the  derivation  of  the  model. 

The  maximum  packing  fraction  of  beads  is  determined  experimentally  as  de¬ 
scribed  in  Ward  et  al.  [WWG09].  The  value  of  </>max  is  measured  to  be  approxi¬ 
mately  0.57-0.58.  Each  preparation  of  particle-laden  fluid  has  a  constant  volume 
of  90  cm3  with  approximately  70  cm3  actually  being  transferred  from  the  jar  into 
the  reservoir.  The  particle-liquid  mixture  is  prepared  according  to  the  desired 
concentration  of  beads  and  silicone  oil. 

To  begin  the  experiment,  the  particle-laden  flow  materials  are  placed  in  a 
plastic  container  and  hand-mixed  using  a  stirring  rod  for  4  minutes,  creating  a 
homogeneous  mixture.  Since  the  density  of  the  particles  is  greater  than  that 
of  the  liquid,  the  particles  settle  out  fairly  quickly  and  the  experiments  must 
be  performed  immediately  after  the  particle-laden  fluid  is  well-mixed.  The  fluid 
is  placed  in  the  reservoir  and  the  gate  is  opened.  A  camera  positioned  above 
the  track  and  perpendicular  to  the  inclined  plane  records  still  images  at  prede¬ 
termined  time  intervals.  The  images  for  the  25  and  30%  particle  concentration 
experiments  are  recorded  at  4  fps  (frames  per  second),  the  35  and  40%  at  2  fps, 
and  the  45%  at  1  fps. 

The  images  are  analyzed  by  an  image-processing  code  and  an  average  front 
position  is  calculated  for  each  image.  Figure  5.2  shows  a  time-series  of  images 
taken  at  45  second  intervals.  In  these  images,  the  development  of  the  fingering 
instability  and  a  dark  particle-rich  ridge  at  the  front  of  the  fluid  can  be  observed. 
Figure  5.3  contains  a  series  of  plots  each  taken  two  minutes  into  the  experiment 
for  a  range  of  particle  concentrations.  Note  that  the  fluids  with  more  particles 
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move  more  slowly,  and  have  a  darker  particle  ridge  at  the  front. 


Figure  5.2:  Time  series  45  seconds  apart,  35%  particle  concentration,  45°  incli¬ 
nation. 


Figure  5.3:  Varying  particle  concentrations  25,  30,  35,  40,  and  45%  at  two  min¬ 
utes. 

The  images  recorded  by  the  camera  are  processed  to  extract  the  profile  of  the 
leading  edge  of  the  fluid.  A  front  position  for  the  central  half  of  the  flow  (away 
from  the  side  walls),  measured  in  pixels  and  averaged  over  approximately  125  data 
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points,  is  calculated  and  later  converted  into  a  physical  distance.  The  average 
front  position  (cm)  is  plotted  against  the  rescaled  time  r  (see  Section  5.2)  to  the 
one-third  power  in  Figure  5.4.  After  an  initial  transient,  the  data  is  approximately 
linear  with  slope  (found  using  least  squares)  analogous  to  the  scaling  constant  C  of 
(5.11),  which  decreases  with  increasing  particle  concentration.  This  qualitatively 
confirms  the  well-mixed  model.  A  quantitative  comparison  of  the  theoretical, 
experimental,  and  numerical  scaling  constants  is  given  in  Figure  5.9.  Note  also 
that  the  data  for  the  different  particle  concentrations  almost  collapse. 


Figure  5.4:  Experimental  results  tracking  the  average  front  position  and  plotting 
it  against  the  rescaled  time  r1/3. 

5.4  Numerical  Simulations 

In  this  section,  we  describe  numerical  simulations  for  the  model  of  a  particle- 
laden  film  with  settling,  i.e.,  (5.7)  and  (5.8).  This  model  requires  a  precursor  as 
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it  has  been  shown  that  solutions  with  settling  depend  singularly  on  the  precursor 
thickness  [CBH08,  ZDB05].  The  simulations  employ  an  upwind  finite- difference 
scheme,  which  is  efficient  for  conservation  laws  with  a  unidirectional  velocity.  We 
use  the  experimental  parameters  described  in  Section  5.3  for  the  viscosity,  density, 
and  the  diameter  of  the  particles.  For  the  maximum  packing  fraction  0max,  we 
use  the  experimentally  determined  value  of  0.58;  see  Ward  et  al.  [WWG09].  We 
assume  initial  data  of  the  step-like  form  described  in  Section  5.2,  with  initial 
height  /3  —  2  and  width  L  =  2.5.  The  precursor  height  for  the  numerical  data 
is  b  =  0.01  unless  otherwise  stated.  The  initial  particle  concentration  is  uniform, 
representing  the  well-mixed  initial  fluid  in  the  physical  experiment. 

We  compute  solutions  with  settling  function  (5.6)  for  m  =  1  and  m  =  5.  Re¬ 
call  that  the  value  for  the  exponent  m  was  determined  experimentally  [RZ54a]. 
For  each  of  the  settling  functions,  Figure  5.5  contains  snapshots  in  time  of  nu¬ 
merical  solutions  of  (5.8).  In  the  left-hand  plot  which  contains  solutions  using 
the  settling  function  (5.6)  with  m  =  1,  the  magnitude  of  the  plots  in  the  vertical 
direction  is  larger  than  solutions  using  the  settling  function  (5.6)  with  m  =  5. 
This  phenomenon  has  been  described  in  Cook  et  al.  [CBH08]  as  a  singular  shock. 
However,  at  these  scales,  the  plots  are  qualitatively  similar.  From  the  modeling 
point  of  view,  the  exponent  m  —  1  seems  to  be  more  appropriate  as  it  shows  a 
stronger  ridge  that  is  also  observed  experimentally.  For  m  =  5,  the  concentra¬ 
tions  vary  less  than  1%.  From  the  data  in  Figure  5.5,  one  can  determine  the 
spatial  variations  in  the  viscosity  by  the  empirically  derived  model  (5.4). 

In  Figure  5.6,  we  compare  the  profile  shapes  for  the  analytical  solution  with 
precursor  (left)  and  the  numerical  solution  with  settling  function  (right).  To 
compare  profile  shapes,  we  rescale  the  data  by  the  maximum  height  in  h  and  by 
front  position  in  x.  We  see  that  the  addition  of  particle  settling  exhibits  depar- 
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Figure  5.5:  Numerical  simulations  of  the  particle  concentration  of  solutions 
of  (5.7)  and  (5.8)  with  initial  concentration  0O  =  0.3,  precursor  height  0.01 
and  for  settling  function  (5.6)  with  m  =  1  (left)  and  m  =  5  (right)  at  times 
t  =  30,  60,  90, 120.  While  the  plots  look  qualitatively  similar,  the  more  singular 
solutions  for  rri  —  1  have  larger  maxima. 


Figure  5.6:  A  comparison  of  scaled  solutions.  On  the  left,  simulations  of 
(5.7)-(5.8)  with  precursor  but  no  particle  settling.  On  the  right,  simulations 
with  precursor  and  settling  with  rn  —  1.  Both  simulations  have  precursor  height 
0.05  and  initial  particle  concentration  cj)0  =  0.4. 
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tures  from  pure  self-similarity,  as  seen  in  the  profile  shapes,  for  times  relevant  to 
us.  Another  comparison  of  the  fronts  is  illustrated  in  Figure  5.7.  It  shows  that 
the  height  of  the  precursor  makes  a  much  more  significant  difference  in  the  front 
speed  than  the  addition  of  particle  settling.  With  a  precursor,  the  front  speed  is 
similar  for  both  the  analytical  solution  and  the  settling  models  of  either  power 
and  faster  than  that  of  the  similarity  solution  with  no  precursor. 


Figure  5.7:  Comparison  of  solutions  to  the  settling  model  (5.7)-(5.8)  with 
lytical  solutions  (5.14)  of  the  Riemann  problem  with  precursor  height  b  = 
and  without  precursor  for  initial  particle  concentration  </>o  =  0.3,  at  time  t 

Figure  5.8  is  the  numerical  version  of  Figure  5.4,  which  is  for  physical  experi¬ 
ments.  It  shows  the  same  near-collapse  of  the  front  positions  in  the  right  scaling. 
The  data  is  ordered  monotone  with  the  particle  concentration.  The  qualitative 
agreement  of  Figures  5.4  and  5.8  is  very  good,  whereas  the  quantitative  agreement 
is  not  yet  fully  reached. 

A  quantitative  comparison  between  analytical,  numerical,  and  experimental 
results  is  provided  in  Figure  5.9.  It  shows  the  scaling  constants  derived  from  the 
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Figure  5.8:  The  front  position  xs(r)  plotted  versus  r1/'3  for  different  particle 
concentrations.  The  data  is  nearly  independent  of  the  choice  of  settling  function. 
It  has  precursor  height  b  =  0.01. 

experiments,  the  numerics,  and  the  theory,  all  divided  by  the  theoretical  con¬ 
stants,  in  dependence  of  the  particle  concentration.  The  solid  line  at  1  represents 
the  theoretical  constants  (5.11)  derived  from  the  analytical  solutions  without 
precursor.  The  theoretical  constants  with  precursor  height  b  =  0.01  derived  in 
Section  5.2  are  bigger  than  the  ones  for  b  =  0.  They  are  similar  to  the  numerical 
data  for  precursor  height  b  =  0.01.  The  numerical  values  for  C  for  the  different 
settling  functions  are  so  similar  that  the  points  are  virtually  indistinguishable  on 
this  plot  and  are  represented  by  a  single  averaged  point.  The  numerical  data  is 
the  same  as  in  Figure  5.8.  The  numerical  (experimental)  constants  are  found  by 
the  slope  of  a  linear  interpolation  of  the  data  in  Figure  5.8  (Figure  5.4).  The  best 
agreement  between  the  numerical  results  from  the  system  (5.7)  and  (5.8)  and  the 
experimental  results  is  for  an  initial  particle  concentration  of  approximately  0.45. 
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We  conjecture  that  this  is  a  balancing  point  between  low  particle  concentrations 
when  settling  to  the  substrate  is  dominant  and  high  particle  concentrations  when 
clumping  and  sliding  behavior  dominates.  We  see  that  numerics  as  well  as  theory 
overestimate  the  scaling  constant  found  in  the  physical  experiment  for  low  con¬ 
centrations.  For  high  concentrations,  the  numerical  and  the  experimental  scaling 
constants  almost  agree. 
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Figure  5.9:  Dependence  of  the  scaling  constant  C  on  the  particle  concentration 
0.  The  theoretical  scaling  constant  is  computed  from  (5.11).  The  numerical  and 
experimental  data  come  from  Figures  5.8  and  5.4. 

5.5  Summary 

We  have  developed  a  new  understanding  of  the  role  of  the  precursor  in  lubrication 
models  for  particle-laden  thin  him  hows,  in  determining  the  position  of  the  front 
of  a  particle-laden  huid.  Settling  of  the  particles  in  the  direction  of  the  how 
primarily  affects  the  prohle  of  the  him  as  a  particle-rich  ridge  develops  at  the 
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front,  and  has  a  less  significant  effect  on  the  speed  of  the  front.  The  scaling  for 
the  front  position  that  is  exact  for  analytical  solutions  with  and  without  precursor 
is  also  appropriate  to  leading  order  for  numerical  simulations  of  the  model  with 
particle  settling,  as  well  as  for  the  physical  experiments. 

The  development  of  a  particle-rich  ridge  at  the  front  of  the  fluid  is  seen  in 
the  numerical  data  and  in  the  physical  experiments,  where  it  appears  as  a  dark 
ridge  at  the  front.  However,  it  is  also  clear  that  the  model  does  not  quantita¬ 
tively  reproduce  this  departure  from  self-similarity  of  the  fluid  profile.  A  separate 
paper,  on  the  transverse  fingering  instability  with  surface  tension  [CAB09]  also 
suggests  that  some  important  physics  is  missing  from  the  lubrication  model,  even 
when  surface-tension  effects  are  included  in  the  mathematics.  One  possibility  of 
additional  physics  to  include  in  this  analysis  is  shear- induced  migration  (e.g., 
[LA86,  LA87b])  which  has  recently  been  shown  to  give  quantitatively  accurate 
predictions  of  phase  transitions  between  settling  to  the  substrate  and  settling  to 
the  contact  line  in  constant-flux  experiments  [C00O8].  It  would  also  be  interesting 
to  incorporate  a  precursor  into  the  physical  experiments. 

The  understanding  of  particle-laden  fluid  flows  and  the  relationship  to  clear 
fluid  and  pure  granular  models  is  still  at  a  preliminary  stage.  Appropriate  mod¬ 
eling  for  low  concentrations,  in  which  settling  to  the  substrate  creates  a  phase 
transition  in  which  part  of  the  fluid  becomes  particle-free  needs  more  modeling. 
The  high  particle  concentrations  that  exhibit  sliding  of  large  clumps  will  require 
yet  another  (if  fluid,  then  non-Newtonian)  model. 
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CHAPTER  6 


Experiments  Compared  with  Equilibrium 

Theory 


6.1  Background 

While  many  studies  address  gravity-driven  clear  liquid  flows  [Hup82,  SD85,  THS89, 
JB92,  KB99,  Kon03b]  and  pure  granular  flows  [NHS89,  SH89,  PDS97],  com¬ 
paratively  fewer  studies  have  centered  on  particle-laden  thin  film  flows  [ZDB05, 
CBH08,  C00O8,  WWG09].  Apart  from  complexities  associated  with  moving  con¬ 
tact  lines,  the  study  of  slurries  also  involves  an  intricate  interplay  between  particle 
settling/migration  and  viscous  fingering  mechanisms. 

The  settling  of  particles  in  quiescent  liquids  and  sedimentation  in  suspensions 
have  also  garnered  significant  attention  [RZ54b,  Bat72,  BM73,  BG082,  DA85, 
SMOO].  For  rigid  spherical  particles,  the  well-known  Stokes’  Law  applies,  neglect¬ 
ing  inertial  effects  of  the  liquid  due  to  the  smallness  of  the  Reynolds  number. 
In  order  to  account  for  the  presence  of  a  large  number  of  identical  particles,  the 
velocity  given  by  Stokes’  Law  is  typically  modified  by  a  purely  empirical  mul¬ 
tiplicative  hindrance  function  which  depends  on  the  particle  concentration,  <fi. 
This  function,  typically  denoted  by  /(0),  has  been  a  matter  of  much  discus¬ 
sion  through  the  decades.  A  so-called  Richardson-Zaki  expression  was  proposed, 
where  /(0)  =  (1— (f>)m,  with  m  ~  5.1,  and  found  to  compare  favorably  with  exper- 
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imental  data  for  moderately  dilute  suspensions  [RZ54b].  For  dilute  dispersions, 
m  (1  —  6.550)  was  suggested  in  Batchelor  [Bat 72].  Other,  more  complex, 
expressions  for  /(0)  have  been  discussed  [BM73,  BG082,  DA85,  SMOO].  In  the 
presence  of  shear,  a  hindrance  function  of  form  /(0)  ~  (1  —  0)  was  shown  to  be 
appropriate  in  Schaflinger  et  al.  [SAZ90]. 

Concentrated  suspensions  of  spherical  particles  have  been  shown  to  behave 
curiously  when  subjected  to  shear.  This  phenomena  was  first  detected  in  ex¬ 
periments  with  a  Couette  viscometer,  where  an  unusual  decrease  in  measured 
viscosity  occurred  during  prolonged  shearing.  The  theoretical  framework  for  this 
phenomena  was  laid  out  by  Leighton  and  Acrivos  [LA87a,  LA87b]  and  subse¬ 
quently  rephrased  in  Phillips  et  al.  [PAB92],  Its  key  element  was  shear-induced 
migration ,  a  diffusive  mechanism  resulting  from  gradients  in  both  the  particle  con¬ 
centration  and  suspension  viscosity,  /i(0).  Net  fluxes  caused  by  these  gradients 
were  deduced  by  considering  irreversible  interactions  between  pairs  of  smooth 
spherical  particles  (for  details,  see  [LA87b]).  In  Phillips  et  al.  [PAB92],  the 
predictions  of  this  model  were  shown  to  be  in  excellent  agreement  with  exper¬ 
imental  data  for  Couette  flows,  and  the  use  of  the  model  was  also  extended  to 
flows  of  concentrated  suspensions  through  cylindrical  tubes.  Recently,  the  model 
was  employed  in  Kim  [KLK08]  to  carry  out  numerical  simulations  for  suspension 
flows  in  more  complex  geometries.  Other  studies,  focusing  on  the  migration  of 
particles  in  pressure-driven  channel  flows  [NB94],  steady  and  unsteady  flows  in 
various  geometries  [MB99],  and  inclined  free-surface  channel  flows  [TM05],  were 
carried  out  using  a  different  approach,  based  on  Stokesian  dynamics. 

More  recent  studies  addressed  particle-laden  thin  film  flows  with  contact  lines. 
Zhou  et  al.  [ZDB05]  reported  on  their  preliminary  experimental  results  for  in¬ 
cline  flows  of  suspensions  of  polydisperse  glass  beads  with  diameter  ~  O(lOOyum), 
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focusing  on  a  single  bead  size.  As  they  varied  the  bulk  particle  concentration, 
0o,  and  the  inclination  angle,  a,  they  identified  three  distinct  settling  regimes: 
for  small  0o  and  a  values,  the  particles  would  settle  out  of  the  flow,  clear  liquid 
film  would  flow  over  the  particulate  bed,  and  the  fingering  instability  resulted; 
for  large  values  of  0o  arid  a,  the  particles  would  move  faster  than  the  liquid, 
leading  to  an  aggregation  of  particles  in  the  contact  line  region,  the  formation  of 
a  particle-rich  ridge,  and  almost  complete  suppression  of  the  fingering  instability; 
finally,  intermediate  values  of  these  parameters  would  lead  to  a  well-mixed  regime 
in  which  a  fingering  instability  occurred.  A  theoretical  model  was  also  derived  in 
Zhou  et  al.  [ZDB05],  based  on  the  Navier-Stokes  equations  for  the  liquid  and  a 
diffusive  model  for  the  particle  concentration,  including  capillary  effects  and  hin¬ 
dered  settling.  A  simplified  version  of  this  model,  which  neglected  higher-order 
capillary  terms  was  studied  in  a  shock  dynamics  framework  [ZDB05,  CBH08]. 
Although  successful  in  explaining  qualitatively  the  formation  of  the  particle-rich 
ridge,  they  did  not  provide  a  quantitative  model  nor  did  they  ever  attempt  to 
model  the  other  regimes.  In  an  effort  to  improve  the  understanding  of  these 
regimes,  Cook  [C00O8]  included  shear-induced  migration  in  his  model.  He  as¬ 
sumed  that  the  outcome  of  particle  settling  is  guided  by  the  balance  between 
shear-induced  migration  and  hindered  settling.  Through  his  steady  state  formu¬ 
lation,  he  derived  a  system  of  ODEs  for  0  and  the  shear  stress,  a.  However,  this 
model  did  not  include  a  significant  hindrance  effect  due  to  the  presence  of  the 
solid  track.  Furthermore,  while  he  found  good  agreement  between  the  model’s 
predictions  regarding  the  well-mixed  regime  and  the  experimental  data,  the  data 
itself  was  from  Zhou  et  al.  [ZDB05]  -  old,  preliminary,  and  rather  limited.  Ad¬ 
ditional  work  in  Ward  et  al.  [WWG09]  focused  on  studying  the  propagation  of 
contact  lines  in  particle-laden  thin  him  hows  experimentally.  Apart  from  varying 
0o  and  a,  particle  size  and  density,  and  liquid  viscosity  were  all  varied  in  order 
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to  examine  their  influence  on  the  front  speed.  It  was  found  that  the  dependence 
of  the  front  position  on  time  was  of  a  power-law  type,  with  exponents  similar  to 
the  one-third  proposed  in  Huppert  [Hup82], 

In  this  chapter,  we  carry  out  a  systematic  study  of  the  settling  regimes  over 
a  range  of  particle  sizes  and  liquid  viscosities.  Through  comparison  between 
our  experimental  results  and  the  predictions  of  equilibrium  theory,  we  uncover 
the  transient  nature  of  the  well-mixed  regime,  where  a  bifurcation  to  either  of 
the  remaining  regimes  eventually  occurs.  In  addition,  our  experimental  results 
clearly  indicate  how  the  particle  size  and  the  liquid  viscosity  affect  the  time  scale 
on  which  we  observe  this  transient  regime. 

In  contrast  to  the  preliminary  experimental  results  from  Zhou  et  al.  [ZDB05], 
we  perform  experiments  with  three  different  particle  sizes  and  two  different  liquid 
types,  and  vary  (f>0  and  a  over  wide  ranges  of  values.  The  liquid  viscosity  and  the 
particle  size  are  found  to  affect  the  width  of  the  region  in  (0O>  «) -space  over  which 
the  transition  between  the  settled  and  ridged  regimes  occurs.  Therefore,  we  show 
that  these  parameters  dictate  both  the  likelihood  of  observing  the  well-mixed 
regime  for  given  0O  and  a  values,  and  the  time  scale  over  which  the  well-mixed 
suspension  is  preserved.  Next,  a  theoretical  model  is  derived.  We  consider  the 
steady  state  of  the  system  where  hindered  settling  balances  the  shear-induced 
migration  of  particles.  Our  modeling  approach  is  similar  to  the  one  in  Cook 
[C00O8],  with  one  important  difference:  we  also  include  the  hinderance  to  settling 
due  to  the  presence  of  the  solid  track.  We  proceed  by  showing  excellent  agreement 
between  the  model’s  predictions  and  our  experimental  results  over  all  ranges  of 
viscosities  and  particle  sizes.  Furthermore,  we  show  how  the  results  of  numerical 
simulations  of  our  model  provide  additional  evidence  for  the  transiency  of  the 
well-mixed  regime. 
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6.2  Experimental  Apparatus  and  Techniques 


Figure  6.1  shows  the  experimental  apparatus  we  use,  which  is  the  same  as  in 
Section  5.3.  A  liquid/particle  mixture  prepared  beforehand  is  poured  into  the 
reservoir  situated  at  the  top  of  the  track  and  the  gate  is  lifted,  allowing  it  to 
flow  down  the  track,  with  the  contact  line  initially  straight.  Here,  we  only  focus 
on  experiments  with  finite,  constant  suspension  volume.  The  evolution  of  the 
flow  is  monitored  using  a  digital  camera,  which  is  positioned  above  the  track  and 
captures  images  of  the  moving  front  at  predetermined  time  intervals,  typically 
0.25-4  seconds.  Using  this  setup,  we  are  able  to  monitor  the  him  motion,  starting 
from  release,  until  the  front  has  reached  approximately  0.6  m  down  the  track. 
Several  fluorescent  lights  are  placed  below  the  track  for  imaging  purposes,  while 
food-coloring  dye  is  employed  to  enhance  contrast.  Images  are  subsequently 
analyzed,  and  each  experimental  run  is  classified,  based  on  the  observed  settling 
regime,  as  either  ‘settled’,  ‘well-mixed’  or  ‘ridged’  (see  Section  6.3  for  details). 


Figure  6.1:  The  experimental  apparatus. 
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Our  experiments  involve  three  different  particle  types  and  two  different  liq¬ 
uids.  The  particles  are  smooth  spherical  glass  beads  (Ceroglass),  and  we  consider 
three  different  diameters:  cl  =  0.143  mm  (‘PI’),  0.337  mm  (‘P2’),  and  0.625  mm 
(‘P3’).  The  standard  deviation  of  the  particle  diameters  is  26%  for  all  particle 
sizes.  For  the  suspending  liquid,  we  use  polydimethylsiloxane  (PDMS)  (Alfa  Ae- 
sar)  in  two  different  kinematic  viscosities:  v  =  10-4m2/s  (‘LI’)  and  10_3m2/s 
(‘L2’).  The  particles  are  heavy,  i.e. ,  pp  >  pi  for  all  particle  and  liquid  types, 
where  pp  and  pi  are  particle  and  liquid  densities  respectively.  Relevant  material 
parameters  are  summarized  in  Table  6.1. 

Suspensions  are  prepared  by  first  weighing  the  particles  and  PDMS  individ¬ 
ually,  pouring  PDMS  into  a  container,  and  then  adding  particles;  slow  manual 
stirring  is  used  until  a  uniform  mixture  is  obtained.  This  procedure  prevents  the 
formation  of  air  bubbles.  Typically,  no  haste  is  required  between  the  preparation 
of  the  suspension  and  its  release  down  the  track  since  uniformity  of  the  mixture 
is  preserved  for  sufficiently  long  time-intervals.  The  bulk  particle  concentration, 
00)  is  defined  as  0 o  =  Ip/R,  where  V  =  Vi  +  Vp  is  the  total  volume  of  the  mixture, 
and  Vi  and  Vp  are  the  liquid  and  particle  volumes  respectively.  Here,  we  focus  on 
V  between  75  ml  and  103  ml. 


z/(m2/s) 

Pi(  kg/m3) 

Pp(kg/m3) 

d(rnrri) 

LI 

o 

1 

966 

- 

- 

L2 

10”3 

971 

- 

- 

PI 

- 

- 

2475 

0.143 

P2 

- 

- 

2475 

0.337 

P3 

- 

- 

2475 

0.625 

Table  6.1:  Physical  properties  of  the  liquids  and  particles  used  in  the  experiments. 
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The  experiments  are  carried  out  in  open  air  and  at  room  temperature  (298  K), 
maintained  by  the  air-conditioning  unit.  The  fluorescent  lights  we  use  for  imaging 
purposes  radiate  heat,  but  the  amount  is  insufficient  to  affect  either  the  viscosity 
of  the  liquid,  flow  dynamics,  or  observed  particle  behavior  in  any  significant 
manner.  The  track,  gate,  and  reservoir  are  cleaned  after  each  experimental  run 
using  a  squeegee  to  remove  the  excess  particulate  and  dust  which  may  accumulate. 
Although  this  cleaning  procedure  does  not  remove  the  PDMS  entirely,  it  ensures 
reproducibility  of  our  experimental  results. 


PI 

P2 

P3 

LI 

Experiment  C 

Experiments  A,C 

- 

L2 

Experiment  B 

Experiments  A,B 

Experiment  B 

Table  6.2:  The  different  liquid/particle  combinations  we  consider.  We  study  the 
manner  in  which  the  viscosity  of  the  suspending  liquid  (Experiment  A)  or  particle 
size  (Experiments  B  and  C)  affects  the  settling  regime. 

We  carry  out  three  different  sets  of  experiments,  conveniently  summarized 
in  Table  6.2.  In  all  experiments,  we  vary  0O  between  0.25  and  0.50,  and  a 
between  20°  and  50°.  In  Experiment  A,  we  consider  medium-sized  particles, 
P2,  and  both  PDMS  types  in  order  to  study  the  influence  of  the  viscosity  of 
the  suspending  liquid  on  the  settling  regime.  Experiments  B  and  C  focus  on 
studying  the  influence  of  particle  size,  by  fixing  the  liquid  type  (L2  in  B  and 
LI  in  C)  and  varying  the  particle  size.  When  LI  PDMS  and  P3  particles  are 
used,  rapid  settling  occurs.  Regardless  of  our  best  efforts,  a  significant  fraction 
of  particles  often  settles  to  the  bottom  of  the  reservoir  before  the  suspension  is 
ever  released  down  the  track.  Hence,  we  omit  experiments  with  this  mixture. 
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6.3  Experimental  Results 


In  all  experiments,  the  observed  flows  are  relatively  slow.  In  addition,  settling 
behavior  can  only  be  classified  after  an  initial  transient  stage  which  typically  lasts 
up  to  900  seconds.  The  settling  regimes  observed  in  our  experiments  resemble  the 
ones  discussed  in  Zhou  et  al.  [ZDB05] :  each  experimental  run  is  labeled  as  either 
settled,  well-mixed,  or  ridged.  Typical  examples  of  these  regimes  are  shown  in 
Figure  6.2.  In  general,  these  three  regimes  occur  in  each  experimental  set,  A,  B, 
and  C  (an  exception  is  discussed  below).  In  the  settled  regime,  the  particles  tend 
to  quickly  settle  out  of  the  flow,  forming  a  particulate  bed,  with  the  suspending 
liquid  moving  down  the  track  faster  than  the  particles.  Virtually  clear  liquid 
him  ultimately  leaves  the  particulate  bed  far  behind  and  develops  the  fingering 
instability  as  described  in  Huppert  [Hup82],  Typically,  this  regime  occurs  for 
small  values  of  0O  and  a.  In  contrast,  when  0O  and  a  are  large,  the  particles 
move  faster  than  the  suspending  liquid  and  aggregate  in  the  contact  line  region, 
forming  a  particle-rich  ridge,  often  several  times  thicker  than  the  trailing  him. 
Hence,  we  refer  to  this  regime  as  ridged.  Large  particle  concentrations  at  the 
front  appear  to  suppress  the  fingering  instability.  Intermediate  values  of  0o  and 
ol  lead  to  the  well-mixed  regime,  where  the  particle  concentration  remains  almost 
uniform  throughout  the  him.  The  fingering  instability  occurs,  but  compared  to 
the  settled  regime,  it  is  typically  characterized  by  a  longer  wavelength. 

For  any  particular  liquid/particle  combination,  we  classify  each  experimental 
run  and  compile  the  results  of  the  classification  in  a  corresponding  (0O,  a)  phase 
diagram  (see  Figures  6.3,  6.4,  and  6.5).  Same-type  settling  behaviors  cluster  in 
these  diagrams,  forming  distinct  bands.  We  use  color  to  label  these  regime  bands: 
white  for  settled,  light  for  well-mixed,  and  dark  for  ridged.  Figures  6.3,  6.4,  and 
6.5  each  include  a  curve  superposed  on  the  experimental  results.  These  curves 
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represent  a  well-mixed  regime  prediction  (i.e.  <j)'  =  0,  where  (j>(z)  is  particle 
concentration)  of  our  theoretical  model.  The  discussion  regarding  the  model  and 
the  agreement  between  its  predictions  and  the  experimental  results  is  given  in 
Sections  6.4  and  6.5. 


a)  b)  c) 


Figure  6.2:  The  settling  regimes:  a)  settled,  b)  well-mixed,  and  c)  ridged.  The 
fingering  instability  typical  of  clear  liquid  flows  is  only  observed  in  the  settled 
and  well-mixed  regimes. 

We  note  that  other  more  complex  behavior  also  occurs.  In  some  experimental 
runs,  we  notice  capillary  motion  of  the  particles  along  the  side  walls  or  their 
alignment  in  linear  streaks  along  the  track.  Irregularities  in  the  shape  and  size 
of  fingers,  and  extreme  cases  of  the  ridged  regime,  where  sections  of  the  contact 
line  experience  jamming  with  particles,  become  solid-like  and  virtually  break 
off  in  blocks  are  also  observed.  These  phenomena  are  attributed  to  either  the 
finite  width  of  the  track  and  the  presence  of  side  walls  or  complex  interplay 
between  particle  migration  and  contact  line  effects.  While  very  intriguing,  we 
leave  detailed  study  of  such  complex  behavior  for  future  work,  and  focus  here 
on  the  three  settling  regimes  described  above.  We  proceed  by  discussing  the 
dependence  of  settling  behavior  on  the  viscosity  of  the  suspending  liquid  and  the 
particle  size  by  presenting  results  of  Experiments  A,  B,  and  C. 
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In  Experiment  A,  we  consider  intermediate  size  particles,  P2,  and  both  the 
low  and  the  high  viscosity  liquid  (LI  and  L2  respectively).  This  allows  us  to 
study  the  dependence  of  observed  settling  behavior  on  PDMS  viscosity.  Since, 
based  on  Stokes’  Law  (e.g.,  [Bat72]),  the  settling  velocity  of  particles  is  inversely 
proportional  to  the  liquid  viscosity,  a  decrease  in  PDMS  viscosity  should  result 
in  an  enhanced  tendency  of  the  particles  to  settle  out  of  the  flow.  Figure  6.3 
shows  phase  diagrams  which  result  for  low  a)  and  high  viscosity  b).  At  first 
glance,  the  outcomes  appear  to  be  rather  similar  for  the  two  liquids,  although 
the  band  for  the  settled  regime  is  somewhat  wider  in  Figure  6.3  a)  compared  to 
b),  confirming  our  expectations  based  on  the  settling  time.  A  closer  inspection 
also  reveals  that  the  well-mixed  band  is  noticeably  wider  when  the  viscosity  of 
the  suspending  liquid  is  lower  (Figure  6.3  a)).  In  order  to  better  understand 
this  difference,  we  consider  the  time  scales  of  the  motion  of  the  front  and  the 
settling  of  the  particles.  For  clear  fluid  inclined  flows,  the  former  time  scale  is 
proportional  to  the  viscosity  of  the  liquid  (e.g.,  [ZDB05]),  T  cc  /i.  However, 
for  suspensions  p  =  p(4>).  For  estimation  purposes,  we  may  assume  uniformity 
of  the  slurry  (i.e.,  well-mixed  case),  and  by  using  p  =  piv{l  —  (j)/(j)m ax)-2  as  in 
[TM05],  we  get  that  T  oc  v.  It  is  now  clear  that  a  decrease  in  the  viscosity  of  the 
suspending  liquid,  with  all  the  other  material  parameters  fixed,  leads  to  a  faster 
propagation  of  the  suspension  front.  On  the  other  hand,  for  purely  gravity-driven 
particle  settling,  the  relevant  time  scale  is  viscous  and  also  directly  proportional 
to  v  (see  Section  6.4)  -  a  decrease  in  viscosity  leads  to  a  faster  settling  of  particles. 
Therefore,  if  gravity  were  the  only  mechanism  responsible  for  particle  settling, 
the  width  of  the  well-mixed  bands  in  Figure  6.3  would  have  been  independent  of 
the  viscosity.  The  fact  it  is  not  suggests  that  some  part  of  the  relevant  settling 
dynamics  occurs  on  a  time  scale  other  than  the  viscous  one.  In  Section  6.4,  we 
conjecture  that  the  settling  behavior  is  governed  by  a  balance  between  settling 
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due  to  gravity  and  shear- induced  migration.  This  balance  may  lead  to  settling 
which  occurs  on  a  different  time  scale,  introducing  a  correction  to  the  purely 
viscous  one.  Accordingly,  while  a  decrease  in  the  viscosity  clearly  affects  the 
motion  of  the  front,  it  may  only  have  a  relatively  minor  effect  on  the  settling 
rate. 


Figure  6.3:  Phase  diagrams  for  Experiment  A.  Particle  type  is  fixed  (P2),  viscos¬ 
ity  of  the  suspending  liquid  is  varied:  a)  low  (LI);  and  b)  high  (L2).  The  symbols 
denote  the  regimes  observed  in  experimental  runs:  circles  (•)  for  settled,  triangles 
(a)  for  well-mixed,  and  diamonds  (♦)  for  ridged.  The  solid  curve  represents  the 
prediction  of  our  theoretical  model  (see  Sections  6.4  and  6.5)  for  a  regime  where 
(f)'  —  0  (well-mixed). 


Next,  we  examine  the  manner  in  which  particle  size  affects  the  settling  regime. 
For  this  purpose,  we  carry  out  Experiments  B  and  C,  where  the  liquid  type  is 
fixed  and  the  particle  size  is  varied. 

In  Experiment  B,  we  consider  the  high  viscosity  suspending  liquid,  L2,  and  all 
three  particle  sizes,  PI,  P2,  and  P3  (note:  experimental  runs  with  the  L2/P2  com¬ 
bination  have  already  been  carried  out  in  Experiment  A).  According  to  Stokes’ 
Law,  the  settling  velocity  is  proportional  to  a2,  and  hence,  the  largest  particles 
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are  most  likely  to  settle  out  of  the  flow.  Based  on  this  reasoning,  the  settled  band 
should  be  widest  for  P3.  The  phase  diagrams  resulting  from  Experiment  B  are 
shown  in  Figure  6.4. 


Figure  6.4:  Phase  diagrams  for  Experiment  B.  Viscosity  of  the  suspending  liquid 
is  fixed  (high  viscosity,  L2),  particle  size  is  varied:  a)  small  (PI);  b)  intermediate 
(P2);  and  c)  large  (P3).  The  symbols  denote  the  regimes  observed  in  experimental 
runs:  circles  (•)  for  settled,  triangles  (a)  for  well-mixed,  and  diamonds  (♦)  for 
ridged.  The  solid  curve  represents  the  prediction  of  our  theoretical  model  (see 
Sections  6.4  and  6.5)  for  a  regime  where  4>'  —  0  (well-mixed). 


Compared  to  Experiment  A,  the  differences  between  the  diagrams  are  much 
more  pronounced.  The  speculation  based  on  Stokes’  Law  is  again  proven  correct 
compared  to  Figures  6.4  a)  and  b),  the  band  corresponding  to  the  settled  regime 
is  widest  in  Figure  6.4  c);  it  is  narrowest  for  the  smallest  particles  in  Figure  6.4 
a).  But,  the  most  striking  feature  here  is  a  complete  absence  of  the  well-mixed 
regime  for  the  largest  particles  in  Figure  6.4  c)  -  all  considered  runs  with  the 
L2/P3  configuration  resulted  in  either  settled  or  ridged  behavior.  In  addition,  we 
notice  that  the  well-mixed  band  is  significantly  wider  for  the  smallest  particles 
in  Figure  6.4  a)  compared  to  the  intermediate  ones  in  Figure  6.4  b).  Hence,  the 
trend  is  obvious:  for  a  fixed  liquid  viscosity,  an  increase  in  particle  size  makes  the 
well-mixed  outcome  less  likely.  This  result  further  supports  our  hypothesis  re- 
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garding  the  transient  nature  of  the  well-mixed  regime.  In  particular,  the  diffusive 
fluxes  of  particles  due  to  hindered  settling  and  shear-induced  migration  are  both 
proportional  to  a 2  (see  Section  6.4).  Since  these  are,  in  our  opinion,  the  two  main 
mechanisms  of  particle  motion  in  the  system  we  study,  the  smallest  particles  PI 
are  moving  on  a  time  scale  much  longer  than  the  larger  particles  P2  and  P3. 
Consequently,  and  as  seen  in  Figure  6.4  a),  many  runs  involving  PI  remain  well- 
mixed  for  the  duration  of  the  experiment.  We  suspect  that  given  a  longer  track 
and  larger  sample  volume,  the  majority  of  these  flows  would  eventually  bifurcate 
to  either  the  settled  or  ridged  regime.  On  the  other  hand,  the  largest  particles 
P3  move  on  shorter  time  scale  compared  to  both  PI  and  P2.  Therefore,  the  flows 
involving  particles  P3  quickly  bifurcate  to  either  settled  or  ridged,  with  the  latter 
regime  likely  for  most  4>o  and  a  values.  The  complete  absence  of  the  well-mixed 
band  in  Figure  6.4  c)  serves  as  an  indicator  of  just  how  rapid  this  process  is. 

Finally,  in  Experiment  C,  we  study  the  influence  of  particle  size  on  the  settling 
behavior  for  low  viscosity  PDMS,  LI.  We  focus  on  small  and  intermediate  size 
particles,  PI  and  P2.  As  discussed  in  Section  6.2,  we  do  not  consider  the  L1/P3 
combination  since  the  particles  in  all  suspensions  of  that  type  undergo  rapid 
settling  while  still  in  the  reservoir.  We  also  note  that  runs  with  the  L1/P2 
combination  have  already  been  carried  out  in  Experiment  A. 

The  results  of  Experiment  C  are  shown  in  Figure  6.5.  The  trend  observed 
in  Experiment  B  is  very  much  noticeable  here  too:  as  the  particle  size  increases, 
the  uniformity  of  the  suspension  is  less  likely  to  be  preserved.  In  particular, 
the  well-mixed  band  is  significantly  wider  in  Figure  6.5  a)  compared  to  the  one 
in  b)  for  larger  particles.  The  explanation  for  this  trend  is  identical  to  one  for 
Experiment  B.  For  small  particles,  the  time  scale  of  particle  settling  is  much 
longer  than  for  larger  ones;  therefore,  uniformity  of  the  suspension  is  likely  to 
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♦o  % 

Figure  6.5:  Phase  diagrams  for  Experiment  C.  Viscosity  of  the  suspending  liquid 
is  fixed  (low  viscosity,  LI),  particle  size  is  varied:  a)  small  (PI)  and  b)  intermedi¬ 
ate  (P2).  The  symbols  denote  the  regimes  observed  in  experimental  runs:  circles 
(•)  for  settled,  triangles  (a)  for  well-mixed,  and  diamonds  (♦)  for  ridged.  The 
solid  curve  represents  the  prediction  of  our  theoretical  model  (see  Sections  6.4 
and  6.5)  for  a  regime  where  <p'  —  0  (well-mixed). 

be  preserved  longer  for  small  particles.  Furthermore,  a  comparison  of  diagrams 
in  Figures  6.3,  6.4,  and  6.5  reveals  that  the  well-mixed  regime  is  more  likely  to 
occur  for  Ll/Pl  than  for  any  other  liquid/particle  combination  we  consider  - 
the  well-mixed  band  in  Figure  6.5  a)  is  by  far  the  widest.  This  is  particularly 
evident  when  Figures  6.4  a)  and  6.5  a)  are  compared  (small  particle  size,  high  and 
low  viscosity  suspending  liquid  respectively).  The  latter  comparison  also  shows 
that  for  the  smallest  particles,  the  influence  of  viscosity  on  both  prolonging  the 
transient  phase  and  making  it  more  likely  for  a  wider  range  of  0O  and  a  values  is 
much  more  pronounced  than  in  Experiment  A. 

To  summarize,  Experiments  A,  B,  and  C  show  that  the  particle  size  partic¬ 
ularly  affects  settling  behavior.  It  dictates  the  likelihood  of  occurrence  for  the 
settled  regime  and  the  time  scale  for  particle-motion  in  general.  The  viscosity 
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of  the  suspending  liquid  also  influences  the  particles’  motion,  and  increasingly 
so  for  smaller  particles.  The  experiments  also  reveal  the  transient  nature  of  the 
well-mixed  regime.  This  is  evident  from  the  manner  in  which  both  the  particle 
size  and  the  viscosity  of  the  suspending  liquid  affect  the  persistence  of  the  well- 
mixed  regime.  We  argue  that  given  a  longer  track  length,  a  majority  of,  if  not  all, 
the  well-mixed  flows  would  bifurcate  to  either  the  settled  or  ridged  regime.  We 
revisit  this  argument  again  when  we  compare  the  predictions  of  our  theoretical 
model  and  the  experimental  results  in  Section  6.5. 

6.4  Theoretical  Model 

We  consider  a  continuum  model  for  the  particle  concentration,  0.  The  dynamics 
of  0  are  described  by  a  conservation  equation  for  particles,  written  in  an  Eulerian 
reference  frame 


—  V  '  (Jbd  T  Jgrav  T  Jcoll  T  Jvisc)  •  (6.1) 

Here  t  denotes  time,  and  D/Dt  =  d/dt  +  v  •  V,  where  v  =  ( u,w ),  u  and 
w  are  the  components  of  the  liquid  velocity  vector  v  in  the  x-direction  (down 
the  track)  and  ^-direction  (normal  to  track)  respectively.  Equation  (6.1)  includes 
hindered  settling  (Jgrav))  and  shear-induced  migration  effects  (Jcou  and  Jvisc).  It 
also  includes  Brownian  diffusive  flux,  Td  =  — HV0.  We  note  that  since  the  Peclet 
number  corresponding  to  our  problem  is  large  (i.e. ,  Pe  =  7a2 / D  ~  O(103),  where 
7  is  the  magnitude  of  the  local  shear  rate),  we  neglect  this  effect.  The  viscosity  of 
the  suspension  is  a  function  of  the  particle  concentration,  /i  =  p,(0).  Here,  we  nse 
the  expression  from  [TM05],  /i(0)  =  /p(  1  —  0/0max)_2,  where  0max  denotes  the 
maximum  packing  fraction,  and  restricts  the  meaningful  interval  of  values  for  0 
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to  [O,0max],  with  the  mixture  becoming  almost  solid-like  as  </>  — *  0max.  Different 
values  of  0max  have  appeared  in  the  literature,  usually  within  the  range  0.57- 
0.68  (e.g.,  see  [C00O8,  WWG09,  BG082,  SM00,  LA87b,  SAZ90,  PAB92,  MB99, 
TM05]).  We  use  the  procedure  described  in  [WWG09]  and  obtain  </>max  =  0.61. 
Also,  m  =  upi. 

The  settling  of  a  particle  due  to  gravity  is  hindered  by  the  presence  of  other 
particles  and  the  solid  track/wall  [ZDB05].  The  net  flux  of  particles  caused  by 
this  effect  is  given  by 


2 a2(f)(pp  -  pi ) 

jgrav  = - X - f(4)w(z)  g.  (6.2) 

y  pi 

Here,  we  use  the  hindrance  function  from  Schaflinger  et  al.  [SAZ90]:  /(</>)  = 
Pi{  1  —  0)/p.(0).  The  presence  of  a  solid  track  at  z  =  0  is  taken  into  account 
through  w(z)  =  A(z/a)2  / \Jl  +  A2(z/aY  [ZDB05];  A  =  1/18  so  that  w(z)  —>  0 
as  z  — >  0,  and  w  «  1  away  from  z  —  0.  We  note  that  this  is  the  main  difference 
between  our  model  and  the  one  derived  in  [C00O8]:  we  include  the  hindrance  w(z) 
in  our  settling  model,  while  in  Cook  [C00O8],  this  effect  is  neglected  altogether. 

The  effect  of  shear-induced  migration  is  included  in  Equation  (6.1)  through 
two  separate  terms,  Jcon  and  Jvisc.  These  terms  are  defined  as  in  Leighton  and 
Acrivos  [LA87b]  and  Phillips  et  al.  [PAB92],  The  net  flux  of  particles  due  to 
irreversibility  of  collisions  between  pairs  of  particles  is  given  by 


Jcoii  =  -A'coii a2  (V>2V 7  +  </>7V0)  .  (6.3) 

On  the  other  hand,  the  net  flux  due  to  gradients  in  viscosity,  p (</>),  is  given  as 


Jvisc  —  Kvisca  <p  7 


1  dp 
p(4>)  d(j) 


V0. 


(6.4) 
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Here,  Kco n  and  KV1SC  are  proportionality  constants  determined  from  experiments. 
We  follow  Phillips  et  al.  [PAB92]  and  use  Kco n  =  0.41  and  A"visc  =  0.62. 

The  fluxes  given  in  Equations  (6.2),  (6.3),  and  (6.4)  are  all  proportional  to  a2, 
a  fact  we  have  employed  in  Section  6.3,  in  our  argument  regarding  the  influence 
of  particle  size  on  the  settling  behavior  of  particles.  It  is  interesting  to  note  that 
in  Equation  (6.3),  the  first  term  suggests  that  even  if  the  particle  distribution  is 
uniform  (i.e.,  V0  =  0),  migration  will  occur  due  to  gradients  in  the  frequency  of 
irreversible  particle  collisions.  This  migration  will  then  induce  gradients  in  0  and 
hence,  the  second  term  in  Equation  (6.3)  is  activated. 

The  governing  equation  given  in  (6.1)  is  accompanied  by  boundary  conditions 
(zero  normal  flux  at  both  z  =  0  and  z  =  h,  where  h  is  film  thickness)  and  coupled 
to  the  Navier-Stokes  equations  for  a  liquid  with  viscosity  /i(0).  However,  in  order 
to  gain  insight  into  the  settling  behavior  of  particles,  it  is  sufficient  to  consider 
Equation  (6.1)  at  steady  state  [C00O8].  Assuming  the  thin  film  is  flat,  the  steady 
state  is  achieved  when  the  fluxes  given  in  Equations  (6.2),  (6.3),  and  (6.4)  balance 
in  ^-direction 


Jgrav  +  Jeon  +  J  vise  0.  (6.5) 

We  also  assume  that  the  flow  is  simple  and  unidirectional,  so  that  7  =  du/dz. 
Henceforth,  instead  of  using  7,  we  revert  to  the  shear  stress,  a  =  /x(0) 7.  By 
scaling  Equation  (6.5)  using  H  as  the  length  scale  in  ^-direction  and  pigH  sin  a 
as  the  scale  for  cr,  and  differentiating,  we  arrive  at  the  following  ODE 


1  + 


2 (/Wise  -  A' 


coin 


K, 


coll 


oc\)  =  —a'<t>  — 


2 pf  cot  a 
9Kcon 


(1  -<j))w(z),  (6.6) 


where  pf  =  (pp  —  pi)/pi,  and  w(z)  is  the  scaled  version  of  w(z).  Since  we  assume 


that  the  him  is  hat,  the  pressure  is  hydrostatic  in  the  suspension  and  the  (scaled) 
gradient  in  shear  stress  is  given  as  [C00O8] 


a'  =  -(1  + 


(6.7) 


By  substituting  the  expression  from  Equation  (6.7)  into  (6.6),  we  obtain 


1  + 


2(A"visc  -  K, 


coll ) 


K, 


coll 


(T(j)'  =  (1  +  Pf(t>)(t> 


2 pf  cot  a 


9  K, 


coll 


(1  —  (fi)w(z).  (6.8) 


Finally,  the  accompanying  boundary  conditions  are  <r(0)  =  (1  +  p/0o)  and  a(l)  = 
0  [C00O8]. 


The  system  of  ODEs  given  by  Equations  (6.7)-(6.8)  and  accompanying  bound¬ 
ary  conditions  may  be  solved  numerically  for  <fi(z )  and  cr(z).  The  numerical  so¬ 
lutions  of  this  system  are  discussed  in  Section  6.5. 


6.5  Predictions  of  Theoretical  Model  vs.  Experimental 
Results 

The  system  of  equations  for  0(z)  and  cr(z)  given  by  (6.7)-(6.8)  is  solved  using  a 
shooting  method.  The  shooting  is  carried  out  from  z  =  0,  with  0(0)  adjusted  in 
order  to  satisfy  a(l)  =  0.  Since  du/dz  =  cr^z)/ p((p(z)),  once  0(z )  and  p(z)  are 
known,  u(z)  is  readily  found,  by  simply  integrating  once  and  using  the  no-slip 
boundary  condition,  w(0)  =  0. 

As  previously  noted  in  Cook  [C00O8] ,  due  to  the  fact  that  a  is  nonnegative,  0 
is  a  monotonic  function  of  z.  This  is  because  crcf)'  in  Equation  (6.8)  is  determined 
by  a  function  of  0  only,  with  a  single  unstable  root  0(a)  in  the  interval  [0,  0max]- 
Hence,  either  0max  >  0(0)  >  0o  and  0(1)  =  0,  or  0(0)  <  0o  and  0(1)  =  0max, 
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corresponding  to  the  experimentally  observed  settled  and  ridged  regimes  respec¬ 
tively.  The  well-mixed  regime  occurs  when  0'  —  0  and  4>(z)  =  0o  for  0  <  z  <  1. 
Setting  0'  =  0,  </>  =  </>o  in  Equation  (6.8)  allows  us  to  obtain  an  expression  for  a 
in  terms  of  0o 


a  =  tan 


-l 


2P/ 


1  - 


(6.9) 


_9Acoii  (1  +  P/0 o)0o_ 

The  curve  corresponding  to  this  expression  is  shown  as  a  solid  line  compared 
with  the  experimental  results  in  the  phase  diagrams  of  Figures  6.3,  6.4,  and  6.5. 
The  prohles  for  cft(z)  and  u(z)  obtained  by  numerically  solving  Equations  (6.7)- 
(6.8)  using  several  representative  combinations  of  0o  and  a  values  are  shown  in 
Figures  6.6,  6.7,  and  6.8. 


Figure  6.6:  Numerical  solution  for  0o  =  0.250  and  a  =  15°:  (a)  particle  concen¬ 
tration,  0(2);  and  (b)  velocity,  u(z).  Note  that  0max  >  0(0)  >  0o  and  0(1)  =  0. 
This  corresponds  to  the  settled  regime. 

Figure  6.6  shows  the  prohles  for  (0O,  a)  =  (0.250,  15°),  corresponding  to 
the  settled  regime  in  all  phase  diagrams  in  Figures  6.3,  6.4,  and  6.5.  From 
Figure  6.6  a),  it  is  evident  this  is  the  scenario  where  0max  >  0(0)  >  0o  and 


0(1)  =  0.  Most  of  the  particles  are  in  z  <  0.5,  after  which  0  decreases  rapidly. 
Effectively,  the  particle-rich  lower  layer  is  covered  by  a  less  viscous  clear  liquid 
layer.  Furthermore,  in  Figure  6.6  b),  the  velocity  increases  sharply  for  z  >  0.5, 
causing  the  clear  liquid  layer  to  flow  faster  than  the  particle-rich  one.  This  is 
equivalent  to  the  regime  seen  in  Figure  6.2  a)  in  which  particles  settle  to  the 
substrate  and  clear  liquid  continues  down  the  track.  In  this  case,  the  prediction 
of  our  model  agrees  well  with  the  experimental  results. 


Figure  6.7:  Numerical  solution  for  0 o  =  0.475  and  a  =  45°:  (a)  particle  concen¬ 
tration,  (p(z)]  and  (b)  velocity,  u(z).  Note  that  0(0)  <  0o  and  0(1)  =  0max.  This 
corresponds  to  the  ridged  regime. 

In  Figure  6.7,  the  profiles  for  0  and  u  resulting  from  (0O,  a)  =  (0.475;  45°) 
are  given.  We  note  that  in  all  our  experiments  with  these  values  of  0O  and  a, 
the  ridged  regime  occurs  (see  Figures  6.3,  6.4  and  6.5).  Figure  6.7  a)  shows,  in 
contrast  to  Figure  6.6  a),  that  0(0)  <  0O  and  0(1)  =  0max.  Therefore,  particles 
aggregate  close  to  the  free  surface  of  the  film,  and  according  to  Figure  6.7  b),  flow 
faster  than  the  more  dilute  lower  layer.  This  behavior  is  typical  in  the  ridged 
regime  observed  in  experiments  and  the  model  predictions  are  again  in  good 


agreement  with  our  experimental  results. 


Figure  6.8:  Numerical  solution  for  0o  =  0.310  and  a  =  45°:  (a)  particle  concen¬ 
tration,  0(z);  and  (b)  velocity,  u(z).  Note  that  0max  >  0(0)  >  0o  and  0(1)  =  0 
still  apply. 

Finally,  for  (0O,  a)  =  (0.310,  45°),  the  resulting  0  and  u  profiles  are  given  in 
Figure  6.8.  With  the  exception  of  the  experiments  with  the  largest  particles  (see 
Figure  6.4  c)),  this  combination  of  0o  and  a  values  leads  to  a  well-mixed  regime. 
From  Figure  6.8  a),  we  see  that  for  most  of  the  film  thickness,  the  particles 
are  uniformly  distributed.  However,  a  close  inspection  reveals  that  this  case  still 
belongs  to  the  category  of  solutions  to  Equations  (6.7)-(6.8)  where  0max  >  0(0)  > 
0o  and  0(1)  =  0,  namely  the  settled  regime.  In  addition,  Figure  6.8  b)  indicates 
that  the  very  thin  layer  of  clear  liquid  at  the  free  surface  still  flows  faster  than 
the  particle-laden  layer  below  it.  While  our  model  is  a  steady  state  one,  one 
could  clearly  see  how,  in  a  dynamic  setting,  the  situation  shown  in  Figure  6.8 
eventually  leads  to  a  settled  regime,  with  the  top  layer  of  clear  liquid  becoming 
ever  thicker  and  flowing  ever  faster  as  the  system  evolves. 

Next,  we  examine  the  agreement  between  the  a(0o)  well-mixed  curve,  given 
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by  Equation  (6.9),  and  the  experimental  results.  We  note  that,  apart  from  the 
inclusion  of  the  hindrance  due  to  the  presence  of  a  solid  track  in  our  model, 
another  important  difference  between  this  study  and  the  one  in  Cook  [C00O8] 
is  that  we  compare  the  predictions  of  our  model  to  the  results  of  much  more 
extensive  experiments.  For  this  purpose,  we  go  back  to  Figures  6.3,  6.4,  and 
6.5.  In  all  diagrams,  except  the  one  in  Figure  6.4  c),  the  curve  lies  completely 
within  the  well-mixed  band,  in  excellent  agreement  with  the  experimental  results. 
For  the  largest  particles  in  Figure  6.4  c),  the  well-mixed  regime  does  not  occur; 
however,  the  curve  overlaps  a  large  section  of  the  border  between  the  settled 
and  ridged  bands  marking  the  transition  between  these  two  regimes.  This  again 
hints  at  the  transiency  of  the  well-mixed  regime.  The  structure  of  Equations 
(6.7)-(6.8)  also  indicates  that  the  well-mixed  regime  is  an  unstable  root  of  the 
system.  Even  if  </>0  and  a  values  are  adjusted  to  lie  exactly  on  the  well-mixed 
curve,  even  the  smallest  of  perturbations  will  eventually  cause  a  bifurcation  to 
either  the  settled  or  ridged  regime.  The  strongest  evidence  for  this  argument  is 
given  by  Figure  6.4  c),  where  the  relevant  settling  time  scales  are  short  enough  so 
that  the  bifurcation  occurs  rapidly  and  the  well-mixed  band  simply  collapses  to 
the  well-mixed  curve.  Other  phase  diagrams  in  Figures  6.3,  6.4,  and  6.5  are  also 
in  line  with  our  argument,  only  the  time  scales  at  which  the  bifurcation  occurs 
are  much  longer  compared  to  the  one  in  Figure  6.4  c),  leading  to  an  observable 
well-mixed  regime.  A  sufficiently  long  experimental  track  would  allow  for  the 
bifurcation  to  occur,  resulting  in  the  eventual  collapse  of  all  well-mixed  bands  in 
Figures  6.3,  6.4,  and  6.5. 

Finally,  it  is  worthwhile  to  emphasize  that  onr  model  is  a  steady  state  one 
while  the  behavior  shown  (e.g.,  in  Figure  6.2  a)),  where  a  clear  him  leaves  a 
particle-rich  sediment  behind  and  develops  fingers,  is  clearly  a  dynamic  process 
which  could  only  be  captured  by  a  more  complete  evolutionary-type  model.  How- 
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ever,  the  settling  behavior,  which  our  simple  model  correctly  predicts,  is  the  most 
crucial  ingredient,  as  it  truly  sets  the  stage  for  these  more  complex  dynamic  pro¬ 
cesses.  Therefore,  our  model  should  be  considered  as  one  of  the  main  components 
of  any  fully  dynamic  model. 

6.6  Summary 

In  this  chapter,  we  focus  on  experiments  with  particle-laden  thin  him  hows  down 
an  incline,  where  the  effects  of  the  viscosity  of  the  suspending  liquid  and  the  parti¬ 
cle  size  are  examined.  We  observe  that  the  settling  behavior  of  particles  proceeds 
in  three  distinct  regimes:  settled,  well-mixed,  and  ridged,  depending  on  the  bulk 
particle  concentration,  0O,  arid  the  inclination  angle,  a.  Our  theoretical  model, 
based  on  equilibrium  theory  where  hindered  settling  balances  shear-induced  mi¬ 
gration,  is  found  to  be  in  excellent  agreement  with  our  experimental  data.  More 
precisely,  its  predictions  for  the  transition  between  the  settled  and  ridged  regimes 
match  the  experimental  observations  exactly  over  all  ranges  of  viscosities  and  par¬ 
ticle  sizes.  Furthermore,  both  our  model  and  our  experimental  results  suggest 
that  the  intermediate  well-mixed  regime  is  a  transient.  In  particular,  our  equilib¬ 
rium  theory  predicts  no  such  regime;  our  experiments  show  how  the  well-mixed 
band  collapses  as  the  relevant  time  scales  are  changed  by  varying  the  viscosity 
and  the  particle  size.  Therefore,  we  argue  that  the  well-mixed  regime  eventually 
leads  to  a  bifurcation  to  either  the  settled  or  ridged  regime. 

Our  experimental  results  indicate  that  particle  size  is  a  significant  parameter. 
The  likelihood  of  observing  the  well-mixed  regime  increases  with  a  decrease  in 
particle  diameter.  The  viscosity  of  the  suspending  liquid  is  found  to  affect  the 
relevant  time  scale  of  the  flow.  For  the  smallest  considered  particle  size,  the  liq¬ 
uid  viscosity  also  significantly  affects  the  likelihood  of  the  well-mixed  regime:  it 
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is  more  prevalent  in  the  case  of  a  less  viscous  suspending  liquid.  A  combination 
of  a  low  viscosity  liquid  and  small  particles  significantly  affects  the  relevant  time 
scales.  The  flowing  him  runs  out  of  track  length  before  any  substantial  distur¬ 
bance  to  the  uniformity  of  the  suspension  is  observed.  We  argue  that  given  a 
sufficiently  long  track,  the  well-mixed  bands  in  phase  diagrams  such  as  those  in 
Figures  6.3,  6.4,  and  6.5  might  eventually  collapse  to  a  well-mixed  line  given  by 
Equation  (6.9),  so  that  only  the  settled  and  ridged  regimes  are  observed. 

The  development  of  a  tractable  theoretical  model  for  particle-laden  thin  films, 
incorporating  all  the  relevant  physical  mechanisms,  is  an  interesting  problem. 
This  model  would  have  to  account  for  the  momentum  conservation  and  continu¬ 
ity  in  the  liquid,  and  include  the  capillary  and  contact  line  effects,  along  with 
those  relevant  to  particle  migration.  Hence,  this  task  is  still  a  formidable  one. 
Instead,  a  simplified  steady  state  model  is  derived  here.  The  model  considers  a 
balance  between  the  hindered  settling  and  shear-induced  migration  of  particles, 
and  it  is  shown  to  provide  useful  new  information  about  the  settling  behavior. 
This  chapter  therefore  makes  a  significant  step  toward  a  fully  quantitative  model 
by  identifying  the  dominant  equilibrium  physics  for  the  flow.  In  addition,  it  also 
implies  further  modifications  required  in  order  to  fully  understand  the  transiency 
of  the  well-mixed  regime  and  the  intricacies  connected  to  the  time  scales  rele¬ 
vant  to  the  front  motion  and  particle  settling.  Our  study  also  raises  interesting 
questions  regarding  the  motion  of  the  contact  lines  and  the  fingering  instability. 
A  more  complete  theoretical  model  would  allow  for  a  comparison  with  the  time 
dependent  experimental  results  from  Ward  et  al.  [WWG09],  regarding  the  front 
motion  in  particle-laden  films.  In  addition,  the  experiments  with  clear  liquid 
flows  in  Jerrett  and  de  Bruyn  [JB92]  showed  that  once  the  fingering  instability 
occurred,  the  exponent  in  the  power  law  (e.g.,  [Hup82])  describing  the  evolution 
of  the  front  position  was  modified.  Carrying  out  a  similar  study  in  the  particle- 
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laden  setting  would  indeed  be  compelling,  especially  since  it  would  also  allow 
for  the  examination  of  the  connection  between  different  settling  regimes  and  the 
wavelength  of  the  fingering  instability. 
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CHAPTER  7 


Dissertation  Summary 


In  Chapter  1,  we  presented  the  motivation  for  the  ADI  scheme  for  the  particle- 
laden  thin  film  flow  equations  as  well  as  outlined  our  work  on  constant-volume 
problems,  equilibrium  theory,  and  experiments. 

Chapter  2  described  the  system  of  equations  for  particle-laden  thin  him  how. 
One  equation  is  related  to  the  him  thickness,  h,  and  the  other  to  the  particle 
concentration,  q b.  The  equations  are  obtained  using  conservation  of  volume,  for 
the  fluid  as  a  whole  and  for  the  particulate  phase.  The  two  velocity  terms,  vav 
and  vrei,  are  the  average  velocity  of  the  liquid  and  particles  and  the  velocity  of 
the  particle  relative  to  the  liquid,  respectively.  A  shear-induced  diffusion  term  is 
added  in  to  correct  for  an  instability  in  the  particle  concentration. 

Chapter  3  gave  the  derivation  for  the  semi-implicit  ADI  scheme  for  the  particle- 
laden  thin  him  how  equations.  ADI  is  used  for  applicable  terms  and  the  remaining 
ones  are  treated  explicitly.  Forward  Euler  is  employed  for  the  explicit  terms  and 
backward  Euler  for  the  implicit  terms.  The  spatial  discretizations  are  done  via 
centered  differencing.  A  moving  reference  frame  is  utilized  to  reduce  the  size  of 
the  domain  needed  for  simulations,  and  the  speed  for  the  frame  can  be  determined 
based  on  shock  theory  for  the  hrst-order  system  of  equations. 

Various  choices  for  the  implementation  of  the  scheme  and  results  from  the 
simulations  were  discussed  in  Chapter  4.  Based  on  simulations,  the  best  choice, 
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in  terms  of  runtime,  is  performing  iterations  at  each  timestep  and  extrapolating 
the  approximate  terms.  The  scheme  parallelizes  well,  with  nearly  eight  times 
speed-up  for  eight  processors.  The  scheme  serves  as  a  means  to  compare  the 
model  to  experiments.  A  mixture  of  silicone  oil  and  glass  beads  is  used  in  an 
experiment,  and  the  results  from  this  are  compared  with  a  numerical  simulation 
with  the  same  parameters.  The  comparison  shows  good  qualitative  agreement, 
but  lacks  some  quantitative  agreement.  We  also  discuss  the  advantages  of  using 
this  scheme  and  future  work. 

In  Chapter  5,  we  compared  results  for  constant- volume  particle-laden  flows 
in  terms  of  theory,  numerics,  and  experiments.  A  similarity  solution  without  the 
settling  of  particles  is  derived  with  and  without  a  precursor.  In  comparison  with 
the  one-third  power  law,  we  observe  qualitative  agreement  from  experiments. 
Testing  a  different  hindered  settling  function,  we  find  that  for  0max  =  0.67,  the 
exponent  rn  —  1,  instead  of  m  =  5,  produces  a  more  pronounced  particle- rich 
ridge,  which  is  consistent  with  the  behavior  seen  in  experiments.  Comparing  the 
effects  of  the  precursor  and  settling,  the  front  speed  of  the  fluid  is  more  influenced 
by  the  precursor  thickness  than  the  addition  of  particle  settling.  For  the  scaling 
constant,  C,  the  best  agreement  between  theory,  numerics,  and  experiments  oc¬ 
curs  at  0o  —  0.45,  which  is  theorized  to  be  the  point  which  separates  the  regimes 
where  the  particles  settle  to  the  substrate  and  where  the  particles  migrate  to  the 
front  of  the  flow,  forming  a  particle-rich  ridge. 

We  examined  the  effects  of  various  liquid  viscosities  and  particle  sizes  on  the 
behavior  of  the  liquid/particle  mixture  in  Chapter  6.  The  size  of  the  particles  has 
an  effect  on  the  settling  behavior,  with  larger  particle  settling  to  the  substrate 
faster.  This  produces  the  appearance  of  large  well-mixed  regions  for  smaller 
particles.  Due  to  the  speed  of  the  settling,  the  time  scale  for  particle  motion 
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is  affected  as  well.  The  viscosity  of  the  liquid  component  influences  particle 
motion,  especially  for  smaller  particles.  The  transient  nature  of  the  well-mixed 
regime  is  observable  as  the  particle  size  increases.  The  equilibrium  theory,  based 
on  the  model  and  numerics,  predicts  that  the  well-mixed  regime  is  an  unstable 
equilibrium  separating  the  settled  and  ridged  regimes. 
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